Current transformer model with hysteresis for improving the protection response in electrical transmission systems

In this paper, a generic enhanced protection current transformer (CT) model with saturation effects and transient behavior is presented. The model is used for the purpose of analysis and design of power system protection algorithms. Three major classes of protection CT have been modeled which all take into account the nonlinear inductance with remanence effects. The transient short-circuit currents in power systems are simulated under CT saturation condition. The response of a common power system protection algorithm with respect to robustness to nominal parameter variations and sensitivity against maloperation is demonstrated by simulation studies.


Introduction
This paper discusses current transformers (CTs) as an example of time-domain simulation of a nonlinear inductor. An overview of the basic principles of a conventional CT is given in section 2. To gain a deeper insight into the nonlinear and time-variant behavior a generalized steady-state analysis is carried out as well. The state-of-theart methods and models of a nonlinear inductive element including hysteresis will be discussed in section 3. Different hysteresis models are referenced focussing on the practicability for a transient CT simulation. Especially the Preisach model of hysteresis is presented and several features are shown. Section 4 covers the CT differential equation in combination with a hysteresis model. Further requirements to the hysteresis model are pointed out. In section 5 the CT model is used to produce test cases for protection device test. The robustness of the particular power system protection algorithms is evaluated under the influence of CT saturation.

Principles of conventional CTs
Current transformers (CTs) are widely used in power plants, substations and in highvoltage as well as in medium-voltage industrial applications. Their purpose is to measure large load currents and also short-circuit currents with well-defined accuracy. Primary currents in the range of a few kiloamperes are transformed into measurable quantities of a few amperes. According to the winding ratio a primary current i 1 is transformed into a secondary current i 2 considering the current error: (1) Figure 1 shows the principle of a transformer. The primary and the secondary circuit are coupled via a ferromagnetic core. In order to achieve a realistic model of the dynamic behavior it is not sufficient to use a static characteristic between the magnetic flux density B and the magnetic field strength H. The utilization of a hysteresis model can yield additional effects like remanence and dynamical phase lag.
It is known that the CT accuracy is influenced by its physical dimensions: • cross-section of the iron core A F e • mean iron core length l F e • width of the gap l Gap The maximum average permeability µ can be accomplished with a closed iron core. Unlike split cores, a closed profile causes a high remanent flux. Commonly used core materials are silicon-based alloys, e.g. F eSi 5 %. CTs with a closed iron core, such as X, P, P X, P S, T P S, T P X can be summarized in a common class referred to as T P X. All T P X CTs operate with very high accuracy within the linear range below the saturation region. The current error ε i as well as the phase displacement δ i [9,10] are small (see figure 4). The dynamic operation at a higher level can drive the transformer into a remanence point. The remanent flux density B R can take values close to 100 % of the saturation flux density B S as can be seen in figure 2. According to equation 2 the dynamical magnetizing inductance L M takes values from close to 0 in the saturation region up to a few Henry in the linear region. The inductance L M can be modeled either from the primary or from the secondary side. For the following considerations the inductance is placed on the secondary side. Commonly there is only one primary winding, so that w 1 equals 1 and the current transormation only depends on w 2 . For reasons of completeness w 1 is respected anyway in the following equations. The quantity of L M considerably depends on the CT class and the magnetic field strength H (operating point) and the corresponding gradient µ = d B /d H.  (i) The magnetizing characteristic will be linearized. The characteristic is linear even for great values of the magnetic field strength H. (ii) The hysteresis loop will be tightened and thus the remanent flux density B R will be decreased.
For this reason the T P Y class of a CT is called linearized core. The maximum remanent flux is limited to 10 %. The T P Z class has a large gap and is therefore called linear core. It shows a negligible remanent flux. The gradient ∆B /∆H which is depending on the operating point in the B-H-plane is represented by dynamical magnetizing inductance in the equivalent circuit in figure 3. The power system can be modeled as an ideal current source due to the fact that there is no feedback from the CT to the power system. The primary losses in terms of copper loss R 1v and leakage field loss L 1σ therefore become invalid. In case of static or digital relays the resistive part of the burden is dominant and the inductive part can be neglected. The internal resistance R CT is commonly given by the manufacturer. It influences the operational overcurrent factor K scco in case of a connected burden Z b ≈ R b if it differs from the nominal burden Z bn . Figure 4 depicts current error and phase displacement. This phasor diagram is only valid for stationary, sinusoidal quantities. Also, it is assumed that the error current I 0 shows these properties. In that case, the error current could be split into an active component and reactive component. However, this assumption cannot hold due to the fact that the magnetizing inductance is nonlinear. In the following, a hysteresis model will be used to simulate the magnetizing inductance in the whole operating range. By exploiting the nonlinear magnetizing current i µ in combination with the dynamic operation of the hysteresis loop, the subdivision into active and reactive loss components is no longer For small time constants the flux in the core can be reduced in short time. The primary current shows an exponentially decaying DC component. This component is damped by the CT, whereby the filter effect is enhanced by a small time constant (particularly T P Z class).

Jiles-Atherton model
The Jiles-Atherton model based on a first-order differential equation is closely related to Curie-Weiss theory of magnetism. It was first introduced in 1983 and it splits the magnetization M into reversible and irreversible portions. it has been proved to work well for stationary input quantities, but shows disadvantages under asymmetrical excitation as for offset minor loops. Even though approaches are available to improve the characteristic [3], still the model is less flexible concerning parametrization.

Chan model
The newer Chan model is used in LTSpice IV ® software, a popular circuit simulator implementing the SPICE kernel. The model computes the Induction B using a simple algebraic equation with only three shape parameters. Asymmetrical input quantities lead to particular instabilities and to non-physical behavior. This disadvantage is fixed in LTSpice IV but the method is patented by the Linear Technology Corporation [5].

Preisach Model
The scalar Preisach model [14] is based on the integration over an area S, which consists of a positive region S + and negative region S − . The Preisach density function p(α, β) must satisfy the condition as well as The input quantity of H(t) leads to a specific proportion between the positive and the negative subarea. The resulting magnetization at each time instance is calculated by the Preisach formalism: At each time instance a magnetization M is determined. This scalar value depends not only on the current model input H itself, but also depends on the magnetic history: The alternating sequence of the magnetization's local extrema is stored in a descending order. Past extrema can be deleted by current values of the magnetic field strength. The memory vector m(h) works like a LIFO-stack * [15]. * Last-In-First-Out where In many cases it is not possible to find an antiderivate due to the fact that the Preisach density function itself has a complex structure [13]. The calculation of the magnetization M with the Everett function is carried out stepwise. The input values for the Everett function are the current field strength h and the last extremum stored in the memory vector m(h): Along the virgin-curve, a reversal point does not yet exist. Therefore the calculation method is applied. The selection of a suitable Everett function (see also [7]) allows the simulation of the three aforementioned CT classes. Both, the remanent flux B R and the coercitive field strength H c as well as the characteristic of the virgin curve are defined by the Everett function.

CT Modeling
The induced voltage on the secondary clamps of the CT and the magnetizing inductance L M (i µ ) are related to the generalized law of induction. Further, Ampère's law in equation 13 can be applied: If the core is closed and no gap exists, the second term in equation 13 will be zero. After summing up all resistors and inductors , and, using the node rule the differential equation of the CT finally takes the forṁ The nonlinear relationship between the magnetic field strength H ∼ i µ and the magnetic flux density B will be determined from the hysteresis model at any time instance. A convenient hysteresis model takes the magnetic field strength H as an input and the flux density B as an output parameter. Internally, the model calculates the magnetization M which is the material-dependent part of the magnetic field strength. It can be written as a relationship between the magnetization M and the magnetic flux density B. As can be seen in figure 7 that there is a feedback loop from the permeability µ to the magnetizing inductance L M . Previously, the permeability µ is calculated as the gradient ∆B /∆H. The model is numerically solved using Matlab/Simulink ® . Consequently, the dynamical system's behavior significantly depends on the nonlinearity and the hysteresis model. For this purpose the Preisach model has been chosen that reflect these physical characteristics. The inclusion of further effects like frequency-dependent eddy currents or excess losses [13] in the CT model is possible in general. Due to the fact that mainly low frequencies in the range of 50 Hz − 60 Hz are investigated these extensions can be neglected.

Computational effort
The Preisach model in its original form is computationally intensive. However, it can be made more efficient using a transform based on Everett's idea [6]. The need to solve the double integral (equation 6) at each time instance is transformed into the solution of a two-dimensional function. Thereby, the computational effort is enormously decreased.

Vanishing initial susceptibility
The susceptibility χ = d M /d H at the origin as well as in reversal points of the magnetic field strength H moves towards zero in the classical Preisach model. This behavior cannot be approved by physical measurements [8]. For this reason, the classical Preisach model is extended by the Product Preisach model [12]. Major hysteresis loop parameters

Total saturation
Within the Preisach model, there is an maximum output limit where S = S + or S = S − according to figure 5. Beyond this certain point, the transducer is completely saturated and no further increase of the flux density B will take place. Depending on the magnitude of the current flowing through the magnetizing inductance this certain point might be passed. The procedure to cause a slightly increase in flux density beyond the saturation point is a multiplication of the calculated magnetization with the weighted Langevin function [1]. This only affects the region of great field strengths and smoothes the transition around the point of saturation.
5. Modeling of power system faults and the impact of CT saturation A power system protection device operates as a real-time signal-processing-unit. It is located in substations and power plants. A protection device measures currents and voltages in the power system and a binary trip signal is given to the circuit breaker if a power system fault is detected in the protected equipment (line, transformer, busbar, etc.). Figure 9 shows the scheme of a protection device and its connections to the instrument transformers, the circuit breaker and an optional communication channel.
In the following, the CT model as described above will be used to process records of power system faults. The digital recording of currents (and voltages) is processed by the CT model and afterwards these signals are applied to two protection devices Siemens 7SD522, which are capable of  The mentioned protection functions require a current signal and must be robust toward CT saturation effects.

Internal faults
An internal fault is located on the protected line (or electrical equipment in general). The protection devices must under any circumstances determine an internal fault even though the current signal is distorted by CT saturation effects. A high voltage power system including the protected line has been simulated ( Figure 10 Depending on the fault distance from each station at least one device has to trip instantaneously. Figure 11 depicts the reaction on a power system fault with and without CT saturation. The L 1 -L 2 fault is located at 20 % of the line in forward direction. The fault inception takes place at an angle of 150°. The CT is simulated with a ratio of 500 A /1 A as a T P X class CT. The initial remanence in both phases is set to −80 %. The two upper binary traces show the reaction of distance protection and HSD [4] without the influence of CT saturation, the current is ideally transformed (dashed curves). The two lower binary traces show the behavior under CT saturation condition (solid curves). It can be seen that the trip command is delayed only by 2 ms in this case.

External faults
An external fault is located outside of the protected line or the electrical equipment in general as can be seen in figure 12. The protection device must not trip under any circumstances even though the current signal is distorted by saturation effects and lead to an apparently increasing differential current. The single phase-to-earth fault shown in figure 13 is measured by both differential protection devices A and B. Station A again is equipped with an ideal CT, whereas in station B a class T P Z CT is used. Because of the huge transformer ratio of 5000 A /1 A the CT does not saturate. It can be seen that within the linear operating region of the CT of class T P Z the DC component is highly damped. As a result the difference current between the two stations (second analog trace) marks a great measurement error. The devices' reactions are depicted in the lower part. The calculated restraint current (green trace) overbalances the increasing differential current (red trace). In addition the fast charge comparison stage I Dif f >> detects the measuring imprecision and an undesired overfunction can be avoided.

Conclusions
The described model is capable of simulating saturation effects in conventional CTs. Already existing real-time power system modeling units commonly integrate a simple B−H characteristic neglecting the hysteresis effect. Due to the low computational effort, the presented model could be implemented in a variety of power system models and test environments. The CT model can process arbitrary primary currents and therefore has the ability to create test scenarios under the influence of CT saturation which are very complex. Especially, scenarios containing several duty cycles can be investigated based on the assumption of a magnetic history or an initial remanent flux. Concepts for the detection of current waveforms affected by CT saturation have to be robust against further effects like transformer inrush, frequency deviation and others. Also the nonlinear inductor model can be ported to further fields of investigation like electrical generators, motors and transformers.