A Methodology for Remote Sensing Inter-Turn Fault Events in Power System Air-Core Reactors, via Simulation of Magneto Quasi-Static Fields in 2D FDTD

We present a numerical methodology to estimate the transient fault currents and to simulate the remote sensing of transient fault information embedded in the magnetic field emissions caused by inter-turn shorts in 60 Hz air-core reactors, thru a magneto quasi-static (MQS) field approximation in the method of Finite-Difference Time-Domain (FDTD) in 2-dimensional (2D) space. The MQS 2D FDTD fields of reactor in normal operation are scaled by correlation against an equivalent circuit model that is derived from application of basic physics principles to parameters of the 3D air-core reactor. The proposed multi-scale quasi-static modeling methodology, based on the reduced c modification, provides fine-feature access down to the single-wire level and can efficiently estimate the transient fault fields and currents due to turn-to-turn short in a reactor with core height in several meters, core diameter in meters, wire diameter in millimeters, and number of turns in the thousands, at 60 Hz; this is accomplished by using computational resources of a typical laptop computer within seconds or minutes, as opposed to days that would be otherwise required without the reduced c modification.


I. INTRODUCTION
The air-core reactor [2] is primarily used as a current limiting device across power networks around the world. It is comprised of coils of wire wound around a hollow cylindrical dielectric (e.g., PVC in the lab, or frame in the field) with air for its core material, as depicted in Fig. 1. To achieve the same level of inductance in lieu of a magnetic core (such as ironcore), an air-core reactor needs many more wire turns than its magnetic-core counterpart. Iron core reactors are filled with dielectric oil, while air-core reactors are seeing increased use in environmentally sensitive areas. Turn-to-turn faults involving a small number of turns are difficult to detect in air-core reactors, due to the limited voltage drop per turn. Methods The associate editor coordinating the review of this manuscript and approving it for publication was Padmanabh Thakur . applied based on voltage or current measurements to detect faults in iron-core reactors have limited sensitivity to fault detection in air-core reactors [3], [4]. A 2D cross-sectional view in Fig. 2 shows the structure during normal operation, where the parameters of a typical air-core reactor in the field can span a range of values; a reasonable estimate is given in Table 1 [5], [6].
If the dielectric insulation of the wires break (e.g., due to overheating, environmental conditions, etc.) and the winding conductors are exposed to form an inter-turn electric short as depicted in Fig. 3, then a small initial fault current can evolve over time to cause more significant damage to the reactor, requiring it to be removed from service.
Previous investigators have proposed models and methods for detecting the fault current due to inter-turn short (also called an inter-winding or turn-to-turn short). For example, VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ FIGURE 1. Conceptual view of a typical 3D reactor with air core (hollow gray cylinder) and winding wire conductors (orange lines), where core radius is r c and core height is h c . This is a conceptual perspective because the parameter scales were chosen to make distinct wires, with d w in millimeters, visible relative to the entire core with r c and h c in meters.
The work by [8] performed simulations of equivalent circuit models based on the finite element method (FEM) to study the frequency-domain variations in spatial distribution of magnetic field. An independent, but somewhat related, area in power networks that often requires electromagnetic modeling and analysis of solenoids is the rotating machines which have been the subject of a lot more analysis [9]- [12] due to their existence in much larger quantities and their higher economic value. Although, reactors and rotating machines are devices with different dimensions and applications, previous works may provide useful insights about the state of computational electromagnetic (CEM) modeling in dense coils of wiring, at 60 Hz. For example, the work by [13] performed fault analysis of motors using a CEM software based on 2D FEM, where each coil has 26 turns and turn-to-turn fault is induced by shorting 6 turns. The work by [14] developed a phase variable model of machines based on FEM combined with wavelet analysis, to predict turn-to-turn short circuit fault, and applied it to a 2-hp/6-pole/36-slot motor where up to 3 turns are shorted.
Our main contribution in this work is to develop an efficient time-domain CEM modeling methodology to sense inter-turn faults based on transient magnetic field emissions, rather than using voltage or current quantities. More specifically, we present a modeling and simulation methodology to remote-sense the transient fault information embedded FIGURE 2. A 2D cross-section view of an air-core reactor. The core height is along the y-axis and the winding wire conductors (orange disks) carry current I z along the z-axis. d c is the diameter of the core, h c is the height of the core, N T is the number of turns, N wl is the number of winding layers, w wl is the total width of the winding layers along the x-axis, r w is the radius of the wire's conductor, r i is the radius of the wire's dielectric insulation (blue disks), d i is the pitch (center-to-center distance) of two nearest wires. The permittivity, permeability, and conductivity of the material is given respectively for the conductor { c , µ c , σ c }, for the insulation dielectric { d , µ d , σ d }, and for the air core { 0 , µ 0 , σ 0 }. Unit vectors in the Cartesian coordinate system are {x,ŷ ,ẑ}. The C E .d l = V is integrated along the contour C v , parallel toŷ at each index m alongx, to compute the voltage V coil across terminals of the coil with N T turns.  Table of parameter values for a typical air-core reactor. Usually in the field for a 3-phase 60 (Hz) power system, three such reactors are connected in a linear or triangular alignment. We define the nominal parameter values as N T = 1000, r c = 1.0 m, I 0 = 1.0 A, d = 0 , σ d = σ 0 , and N wl = 1; unless otherwise explicitly stated, the nominal values are assumed throughout this paper.

FIGURE 3.
Fault current can be induced if wire conductors (orange disks) are shorted (red rectangle) across layers or within a layer. This can occur by wear-out in wire insulation due to repeated surges, environmental effects, etc. Ampere's law C H · d l = I is integrated along the closed contour C s to compute the fault current on the short conductor.
in the magnetic field emissions, due to inter-winding shorts in air-core reactors with realistic parameters specified in Table 1, where core height is in several meters, core radius is in meters, wire diameter is in millimeters, and number of wire turns is in the thousands. The estimate is computed within seconds or minutes, using computational resources of only a typical laptop computer. The time-domain CEM model is implemented in the method of finite-difference timedomain (FDTD) in the magneto quasi-static (MQS) regime, with fine spatial resolution that provides access to the smallest feature (single-wire) all the way to the largest feature (core height) in one model. We believe the proposed methodology for modeling transient fault events in air-core reactors is unique and has not been presented in the literature previously, to the best of our knowledge.
We develop the proposed modeling methodology in the remainder of this paper which is organized as follows.
In Section II, we provide an overview of the essential challenges in using the original FDTD algorithm to solve the electromagnetic (EM) fields of a multi-scale and quasi-static problem. In Section III-A, we rely on basic physics principles to develop an approximate equivalent circuit model of the 3D reactor structure. In Section III-B, we define the excitation current source as a Gaussian pulse with a 60 Hz half-bandwidth at −3 dB, and derive the time-domain voltage response of the equivalent circuit model for later use in correlation and scaling of the FDTD model. In Section III-C, we apply the reduced c method, through modification of free-space permittivity 0 , to the standard Yee FDTD algorithm and obtain an MQS 2D FDTD model to estimate the electric field E(t, x, y) and magnetic field H (t, x, y) in 2D space (in the x-y plane) as a function of time t (s). In Section III-D, we use the voltage and current approximations from the circuit model in normal operation to scale and correlate the mathematical EM fields from the MQS 2D FDTD model. In Section IV-A, we correlate the FDTD model against the equivalent circuit model via several numerical experiments. In Section IV-B, we simulate the transient fault currents on inter-winding shorts, and in Section IV-C, we demonstrate the proposed modeling methodology by computing the time-varying magnetic field emissions to simulate remote sensing a transient fault event in the reactor of Table 1. We conclude with closing remarks in Section V.

II. FDTD FOR MULTI-SCALE & QUASI-STATIC REGIME
Given that 60 Hz is a relatively low frequency, it may be tempting to pursue a static solution; however, such a solution is inadequate for quantifying the transient nature of the fault currents and fields. The structure depicted in Fig. 2 with parameter values specified in Table 1 does not lend itself easily to Yee's original FDTD method [15]- [17] for solving the transient EM fields, due to its multi-scale and quasi-static nature; as explained below.
There is significant spatial disparity between the structure's size and the operational wavelength. Specifically, the smallest feature of the structure (d w = 4.12 mm) is about nine orders-of-magnitude smaller than the operational wavelength (λ 0 = 5000 km). In addition, the smallest feature is about three orders-of-magnitude smaller than the largest feature (h c 3 − 8 m) that further increases the computational burden when using a uniform grid, since the bottleneck on the spatial discretization of the largest feature (h c ) is dictated by the smallest feature (d w ).
In pursuit of a transient solution, starting with the standard Yee's 2D FDTD method, a typical spatial discretization based purely on wavelength, may be x = λ 0 /20 = 250 km [18]; however, our goal of modeling inter-turn shorts requires access to the single-wire level, and a x in kilometers clearly fails to provide sufficient spatial resolution for sampling each wire in millimeters. A more reasonable choice for x may be some fraction of the wire diameter; e.g., x = d w /4 = 1.03 mm; however, considering the Courant stability criteria (i.e., t < .0E+8 m/s), that choice of x dictates an extremely small t = 2.43 ps. Considering that at 60 Hz, one period T 0 = 1/f 0 16.667 ms, it would require approximately 6.87E+9 time-steps to simulate the response to just one period T 0 of a time-harmonic excitation; that is an enormous number of time iterations which exerts increasingly exorbitant computational costs across 1D, 2D, and 3D FDTD on any ordinary, yet respectable, workstation computer today.
The issue of massive time-iterations in the quasi-static regime, due to tiny temporal discretization t in the Courant-constrained central-difference based FDTD, has been studied by previous investigators [19]- [24] for various other applications. The basic idea is that at low frequencies (e.g., 60 Hz) where the wavelength is extremely large compared to the (electrically small) structure's size, the wave's propagation-time across the characteristic length of the structure is much smaller than the time period of interest T 0 (i.e., h c /c 0 T 0 ) [25]; yielding conditions that are suitable for quasi-static approximation of Maxwell's (∇ × E and ∇ × H ) equations which, under the above circumstances, are weakly coupled, or almost (but not completely) decoupled.
An interesting idea, called reduced c, was proposed by previous investigators [19]- [21] for various different applications, where the free-space material properties are modified to reduce the propagation velocity and thereby increase t closer in order-of-magnitude to T 0 . For this paper, we utilize the reduced c method, due to its relative simplicity and convenience in requiring minimal modification to the original Yee's FDTD algorithm.
The 2D FDTD simulation results presented in section III, section IV-A, and section IV-B were generated using a fairly typical (by today's standards) laptop computer running 64bit Linux [26] operating system on a dual-core CPU [27] and random access memory (RAM) of 24.0 GB. Each 2D FDTD simulation that was run at the fine (coarse) sampling corner in the range of spatial (8) and temporal (9) discretization, generally took minutes (seconds) to estimate the EM fields of the nominal 2D air-core reactor in Table 1; in contrast, without the reduced c modification each FDTD simulation run would have required days to complete.

A. EQUIVALENT CIRCUIT MODEL
The circuit model is derived from basic physics principles and the 3D parameters of air-core reactor in Fig. 1 and Table 1.
Each loop of the coil may be modeled as an equivalent circuit comprised of a resistor in series with an inductor [28], where R p m ( ) is the resistance per loop (1), and L p m (H) is the inductance per loop (2); the winding layer number is designated by the subscript index m ∈ {1, 2, 3, . . . }. Internal inductance of the wire is considered negligible here, and thus ignored.
where the resistance per-unit-length (pul) of the wire is defined as Subsequently, the impedance Z coil ( ) across terminals of the entire coil in normal operation, with N T turns and N wl winding layers, may be obtained by where the imaginary number  = √ −1, the angular frequency ω = 2πf (rad/s), and the cyclic frequency is f (Hz).

B. CURRENT SOURCE EXCITATION, AND VOLTAGE RESPONSE
The reactor model is excited by an ideal current source that has a Gaussian pulse waveform with a 60 Hz half-bandwidth at −3 dB, and a peak value of I 0 (A), given by where t pk is the time at which the pulse peak occurs, and t sp is the pulse time-spread. To achieve a 60 Hz half-bandwidth at −3 dB from the peak, we set t sp = 2.20 ms while ensuring t pk is greater than approximately 5 t sp ; i.e., t pk 11.0 ms. Given that in our case, skin-effect [29] is negligible (i.e., at 60 Hz, copper has a skin depth δ s √ 2/(ωµσ ) = 8.53 mm > d w ), the source current may be approximated to have a uniform distribution across each wire's conductor region.
Multiplying the Fourier transform of (4) by (3), then taking the inverse Fourier transform of the product [30], yields the time-domain voltage response (5) of the entire coil to the Gaussian pulse excitation in normal operation.
Assuming an infinitely long wire in theẑ direction and a source current having only a z-component I z (ρ) flowing alongẑ, permits a 2D approximation where it is convenient to assume z-invariance (i.e., ∂ ∂ z = 0) to formulate a transversemagnetic-to-z (TM z ) problem [29] with the z-component of H-field H z = 0; this approximation eliminates three field components E x , E y , H z and one spatial variable z, and leads to the modified Maxwell's equations (6) with only three non-zero field components E z , H x , H y , where the electric field has only a z-component E z alongẑ, and the magnetic field possess only an angular component H φ alongφ [29]. where the electric flux density D = 0 r E (coulombs/m 2 ), r is the relative permittivity of the medium, and the magnetic flux density B = µ 0 H (webers/m 2 ). For the purposes of this work, an electromagnetic modeling approach in 2D space was presumed sufficient, and preferred over a 3D approach which would substantially increase computational costs. The 2D environment assumes that the inter-winding short conductor is a rectangular slab laid in the x-y plane that extends infinitely alongẑ, and induces a transient fault current alongẑ that is attributed to the magnetic field components H x , H y (or H φ alongφ).
The MQS approximation (7) is implemented in 2D FDTD. Furthermore, it was found that uniform spatial discretization in the range of (8) provides sufficient spatial resolution. Considering the Courant stability criteria in 2D FDTD, that choice of x leads to uniform temporal discretization in the range of (9) across each time-step n.
t ∈ [0.0605 ms, 0.121 ms] Given the goal of this work is to present a CEM modeling methodology for numerical estimation of the transient fault currents and remote-sensing transient fault events due to inter-winding shorts in air-core reactors, further details of formulation and implementation of MQS equations in 2D FDTD fall outside the scope of the present paper; however, those details may be submitted as a separate publication in the future.  Fig. 5, we compare the voltage and current from MQS 2D FDTD vs. the equivalent circuit model. As can be observed, the normalized (relative) version of voltage waveforms in Fig. 5(a),(b) are in good agreement; a similar assessment is made of the normalized (relative) version of current waveforms in Fig. 5(c),(d). However, in absolute terms, there is about five orders of magnitude discrepancy in voltage amplitude values between the un-scaled MQS 2D FDTD solution vs. the equivalent circuit solution. The amplitude discrepancy across scaled vs. un-scaled plots of the current appear to be less significant than that of the voltage; however, the current is computed for one single wire turn of the coil, while the voltage is computed across all N T wire turns of the coil.  Table 1.

1) SCALING CURRENT FROM MQS 2D FDTD
Ampere's law [18], [29] may be applied to the H-field computed in FDTD to find the current I z through the line VOLUME 8, 2020 integral of H along a closed contour C i around a single wire; i.e., C H · d l = I z , see Fig. 4. The current I z3 is computed along contour of radius 3 × r w from center of a single wire. The MQS 2D FDTD current solution I FDTD may be scaled to the equivalent circuit solution I FDTD scaled , by using the current scaling factor S I in where S I is a constant based on assumptions in the MQS 2D FDTD; e.g., in this case for nominal parameter values of Table 1, it is empirically determined that S I 2.0. Scaled current results are shown in Fig. 6, where I s was aligned with the curve of r w = 0.81 mm for reference; notice a slight time-delay that increases with increasing r w , given that the integration contour radius increases with r w . Additional small discrepancies may be due to approximating a circular contour in a rectangular FDTD grid with finite spatial discretization.

2) SCALING VOLTAGE FROM MQS 2D FDTD
Maxwell-Faraday's law [18], [29] describes the relation between the E-field obtained from FDTD and the voltage V yz through the line integral of E on the contour C across terminals of the coil alongŷ over N T turns; i.e., C E.
where magnetic flux density B = µ 0 H , and V is the voltage around the closed contour C.
Referring to Fig. 2, the voltage V coil across terminals of the coil with N T turns, is computed by C E · d l = V coil along the contour C v , according to (11). The electric field E m , at each winding layer location index m, is integrated in the dielectric region y ∈ [y j_s , y j_e ] between two conductors, and summed over N T wire turns and N wl winding layers.
where j is the location index of position y, and integration for each wire turn begins at y j_s and ends at y j_e . The MQS 2D FDTD voltage solution V FDTD may be scaled to the equivalent circuit solution V FDTD scaled by using the voltage scaling factor S V in (12), with S V given in (13), wherein the assumption d c r w is implicit.
where d w = d w / x, and functions α( d w ) and κ(r w ) are coefficients dependent on assumptions in the MQS 2D FDTD, and empirically determined through numerical experiments. Fixing the reactor at the nominal values in Table 1 Table 1. Fig. 7, where the FDTD solution slightly overshoots the equivalent circuit solution; this may be explained by a lower overall inductance and resistance implicit in the 2D FDTD assumptions, compared to the equivalent circuit solution which was derived based on 3D reactor parameters.

The scaled FDTD voltage result is shown in
Note that the current scaling factor S I in (10) is a constant, and for fixed values of d w , r w , d c the voltage scaling factor (13) is also a constant. We assume that the scaling factors are independent of the operation mode, whether normal or faulty; thus, we determine the scaling factors under normal reactor operation and use them to predict the fault currents under inter-winding short conditions. Although a relatively simple equivalent circuit model, based on 3D reactor parameters, was used here to scale the mathematical EM fields from the MQS 2D FDTD model, one may instead use measurements of an actual air-core reactor, in the laboratory or on the field, to scale the MQS 2D FDTD model.

A. CORRELATION OF FDTD MODEL AGAINST EQUIVALENT CIRCUIT MODEL, THRU NUMERICAL EXPERIMENTS IN SIMULATION
With the MQS 2D FDTD simulation model scaled to the equivalent circuit solution, we perform five numerical experiments in which we plot V coil (t) vs. time t across variations in r w , N T , r c , N wl , and d to validate the MQS 2D FDTD model by correlation against the equivalent circuit model.
In the 1 st numerical experiment, we vary the radius of core and plot the voltage while fixing the remaining parameters at nominal values in Table 1; results in Fig. 8 show very good agreement between the FDTD and the equivalent circuit models. In the 2 nd numerical experiment, we vary the number of turns and plot the voltage while fixing the remaining parameters at nominal values in Table 1; results in Fig. 9 show quite good agreement between the FDTD and the equivalent circuit models.
In the 3 rd numerical experiment, we vary the wire radius and plot the voltage while fixing the remaining parameters at nominal values in Table 1; results in Fig. 10 show good overall agreement between the FDTD and the equivalent circuit models.
In the 4 th numerical experiment, we set N wl = 3, vary N T , and plot the voltage while fixing the remaining reactor parameters at nominal values in Table 1; results in Fig. 11 show reasonable correlation in trends between FDTD and equivalent circuit solutions; however, generally the FDTD  solution peaks at a higher value and decays faster, compared to its equivalent circuit counterpart.
It's worth noting that any discrepancy in results between the FDTD and the equivalent circuit solutions is likely due to the various simplifying assumptions stated previously, including our rather straightforward approach to developing an equivalent circuit model which was simple, yet sufficient for purposes of this work; however, if required, a 2D FDTD solution may be used with a cellular approach [31] to synthesize a more accurate equivalent circuit model which incorporates both the temporal and spatial characteristics of the reactor more rigorously. Furthermore, a 3D FDTD solution may be attempted if additional accuracy is required; however, at substantially increased computational cost.
In the 5 th numerical experiment, we set d and σ d to PVC, vary r w , and plot the current while fixing the remaining VOLUME 8, 2020 FIGURE 11. Voltage V coil (t ) from equivalent circuit (lines) and from FDTD (shapes). With N wl = 3, we vary N T across 100, 500, 1000, for a total number of turns equal to 300, 1500, 3000, respectively.

FIGURE 12.
Wire insulation material is set to PVC. Current from equivalent circuit (solid black), and from scaled FDTD (shapes). The FDTD current is computed via the Ampere's law integration of the magnetic field ( C H.d l = I z ). The source current is I s .
reactor parameters at nominal values in Table 1; results in Fig. 12 indicate a relatively larger time-delay across increasing r w compared to Fig. 6, as changing the wire insulation from air (with r = 1.0) to PVC (with r = 4.0) reduces the wave propagation velocity and increases the capacitance in the dielectric region surrounding each conductor, thereby affecting the temporal and spatial distributions of fields at each time-step n.  {r c + r i , −h s /2} and {r c + r i + 3 r w , +h s /2}. Current is computed along a fixed rectangular contour region C s to the right of the wire conductors (see Fig. 3), of width 2 r w and height 100 d i .

FIGURE 14. The short is a rectangular region with corners at
{r c + r i , −h s /2} and {r c + r i +4 r w , +h s /2}. Current is computed along a fixed rectangular contour region C s to the right of the wire conductors (see Fig. 3), of width 3 r w and height 100 d i .
The short current is computed by application of Ampere's law integration C H .d l = I z along the closed contour C s around the short rectangular conductor. To enable a logical comparison, the current is computed over a fixed area across the five cases of N s , where the difference between normal vs. faulty current is I z = I z | Faulty −I z | Normal ; results are shown in Fig. 13 and Fig. 14 for nominal parameters in Table 1. As can be observed, increasing h s increases I z ; such deviations in the current profile may be exploited to detect a potential fault event due to inter-turn shorts. It is worth noting that direct measurement of the current on the shorted conductor region may not be practical in actual reactor hardware, as the process may be fraught with various technical and safety challenges; instead, we propose remote sensing the fault event through the magnetic field emissions, as described in Section IV-C.

C. SIMULATION OF REMOTE SENSING A TRANSIENT FAULT EVENT, THRU COMPUTATION OF THE TIME-VARYING MAGNETIC FIELD EMISSIONS
As a demonstration of the proposed CEM modeling methodology, we compute the time-varying magnetic field emissions to simulate remote sensing a transient fault event caused by inter-turn shorts in an air-core reactor at nominal parameter values given in Table 1. Two sets of MQS 2D FDTD simulations are run, as the time-varying magnetic field is recorded at several points in vicinity of the reactor in Fig. 15; in set 1 we simulate the reactor under normal operating condition where N s = 0, and in set 2 we simulate the reactor under transient fault event for five short conditions N s = {1, 5, 10, 50, 100}. For each simulation, we record the magnetic field components {H x , H y } vs. time-step, at the red points depicted in Fig. 15. A similar remote-sensing procedure may be performed by placing magnetic-field sensors near an actual reactor in the laboratory or on the field.
The results are shown in Fig. 16 thru Fig. 25. Each figure shows the H-field vs. time in solid black for normal operation. The difference between the H-field in normal operation vs. faulty operation (16) is displayed in colored dotted/dashed lines.
where u ∈ {x, y}.  Several observations are made: (1) as anticipated, clearly there are variations in H-field signal pattern across normal vs. faulty operations, (2) in most cases the faulty signal deviates notably from the normal signal at the same time-step and location in space, where the deviations generally increase with N s , (3) in some cases the H-field signal deviations exhibit symmetric behavior across symmetry lines of the reactor; e.g., the sign of H x flips across the right bottomto-top (Fig. 17, 18, 19) and across the far right bottom-to-top (Fig. 20, 21, 22), while the sign of H y flips across left-toright (Fig. 16, 21), and (4) in other cases the H-field signal deviations exhibit asymmetric pattern relative to location of short; e.g., across the top left-to-right (Fig. 23, 24, 25).
Finally, we set N T = 1000 and N wl = 3 to achieve 3000 total turns, and insert inter-turn short between winding layers m = 3 and m = 2, as shown in Fig. 26. We repeat the above numerical experiments and record the H-field at the      ambient magnetic field; this result makes sense intuitively, given that we increased the number of healthy turns compared to shorted turns.   In light of the above observations, the time-varying H-field emissions across normal vs. faulty operations may be exploited to remotely sense and detect transient fault   events due to inter-turn shorts. The simulation methodology described herein may be used to quantify the magnetic field profile as a function of space, time, and reactor   parameters, to inform the placement and sensitivity of sensors. A possible detection strategy may include exploration of signal processing and machine learning algorithms

V. CONCLUSION
We presented a CEM modeling methodology for estimation of transient fault currents and remote sensing of fault events due to inter-turn shorts in 3D air-core reactors, through a magneto quasi-static approximation in 2D FDTD. An equivalent circuit model of the 3D reactor in normal operation was developed to correlate and scale the MQS 2D FDTD model. The FDTD model was used to compute the transient fault current due to inter-turn shorts. As a demonstration of the proposed CEM modeling methodology, we simulated remote sensing a transient fault event by computing the H-field emissions of a reactor with realistic parameters, compared the transient fault characteristics embedded in the magnetic field signals at several locations around the reactor over time, and suggested some possible strategies for detecting fault events based on signal deviations across normal and faulty operations.
Some topics of future research may include (1) applying signal processing and machine learning algorithms to the emitted H-field signals for real-time detection of transient fault events and locations in the reactor, (2) extending the proposed methodology to MQS 3D FDTD, (3) exploring methods other than FDTD; e.g., time-domain finite element method (FEM) or transmission line modeling (TLM) method, and (4) modeling magnetic-core reactors with hysteresis, where µ r is a non-linear function of H .

ACKNOWLEDGMENT
The 2D FDTD computer program was originally developed for modeling optical interconnects at 200 THz, under the NSF research grant [32]; later, the code was applied to solve the MQS fields of the air-core reactor at 60 Hz, as described in this paper. The authors would like to thank Dr. Law for insightful technical discussion on the merits of modeling inter-turn shorts in 3D reactors thru a 2D approach (J. Law and A. Zadehgol, Personal Communication, March 18, 2020), and Dr. Juan Marulanda Arias for his efforts to research into potential existing alternative solutions, including efforts to model and simulate the structure in ONELAB [33] based on frequency-domain FEM, and discussions around potential application of machine learning to fault detection (J. Marulanda Arias, Personal Communications, July 29, 2019 through May 13, 2020).