Numerical Inductance Calculations Based on First Principles

A method of calculating inductances based on first principles is presented, which has the advantage over the more popular simulators in that fundamental formulas are explicitly used so that a deeper understanding of the inductance calculation is obtained with no need for explicit discretization of the inductor. It also has the advantage over the traditional method of formulas or table lookups in that it can be used for a wider range of configurations. It relies on the use of fast computers with a sophisticated mathematical computing language such as Mathematica to perform the required integration numerically so that the researcher can focus on the physics of the inductance calculation and not on the numerical integration.


Introduction
Inductors are used to provide filtering or energy storage within many types of electrical systems. They are often coiled conductors whose time changing currents induce voltages either in the conductor itself or nearby conductors. They are used extensively in analog circuitry and are an essential component of every electric utility power grid. For most of the twentieth century, calculating inductances was done mainly through the use of formulas or lookup tables [1,2] because the inductance calculation which contains an integrable singularity had been intractable. For the past thirty years or so, simulation tools such as FastHenry [3,4], and Ansys Maxwell [5], have been used to calculate inductances. However, while FastHenty is open-sourced, the cost of Ansys Maxwell or other tools may be beyond the user's means. Also, for these tools, the user is unlikely to be familiar with their algorithms, so the researcher is less likely to develop an understanding of the fundamentals of inductance calculations; moreover an additional discretization of the inductor may be needed. Now, with the use of modern high speed computers that use sophisticated mathematical computing languages, the numerical integration is tractable and so inductances can now be computed using first principles, using a multi-functional computational program that may already be available to the user. With this method, the workings of the numerical integration are performed by the computational program, and so there is no need for the user to become an expert in sophisticated numerical integration techniques or in methods for discretizing the geometry. This method is particularly useful for the design of air core reactors, which protect equipment from potentially damaging power transients. Air reactors, which can consist of tens or even hundreds of interlaced coils, are often custom-designed for a particular operating environment, and inductance simulations of these reactors can be challenging for simulators that were designed with VLSI integrated circuits in mind. Therefore, we studied a method that uses the computational program Mathematica [6] to determine inductances.

Materials and Methods
A formula for the self-inductance of a conductor is derived by computing the magnetic stored energy using the magnetic vector potential [7] which yields, where m o~4 p10 {7 H/m; the primed variables represent the current source points and the unprimed variables, the current field points; the integral is performed over both sets of variables throughout the volume of the conductor; r is the distance between volume elements and f 1 represents the ratio of the current density J 1 to the current I 1 , f 1~J1 =I 1 . The mutual inductance M between two conductors has a similar formula except that the field points are in one conductor and the source points in another, Although there are few highly symmetric configurations that allow Eqns. (1) and (2) to be integrated analytically, these equations can be integrated numerically, in theory, for any configuration of conductors, with the advent of modern computers and techniques.

FastHenry, Ansys Maxwell, and Mathematica
FastHenry determines inductances (and resistances) by approximating a conductor as a a series of discretized rectangular filaments, each having a lumped resistance and inductance. Using mesh analysis, it derives a set of complex linear equations for the filaments (whose number may be in the thousands) and uses advanced numerical tools such as GEMRES [8], along with a multipole approach [9] to approximate the integrals of (1) and (2) that are performed for each filament. The multipole approach uses multipole expansions for (1) and (2) which may be valid for large r and can give an accurate approximation of the integral with far fewer computations. The user specifies the discretization by inputting the coordinates and dimensions of the rectangular filaments.
Ansys Maxwell [5] does not solve Eqns. (1) and (2) but solves Maxwell's differential equations with boundary conditions and user-specified initial conditions using a finite element method in which the solution is numerically obtained for an arbitrary geometry by breaking it down into tetrahedrons, and using second order quadratic expressions as basis functions to represent the electric or magnetic fields, at the vertices and the midpoints of selected edges. The coefficients of the basis functions are solved for using standard matrix techniques. Since Ansys Maxwell is a large finite element program for general electromagnetic problems, the overhead to use it in terms of cost and training is high.
We use Mathematica to determine the inductance by having Mathematica numerically integrate (1) or (2) over the entire volume of the conductor. Mathematica uses a global adaptive strategy, that tries to reach the required precision and accuracy goals of the integral estimate by recursive bisection of the subregion with the largest error estimate into two halves, and computes integral and error estimates for each half. The user inputs a Mathematica file which describes the integral, and does not need to perform any discretization. Mathematica determines the number of recursive bisections to determine an accurate result.
To demonstrate the accuracy of this technique, we use configurations with known inductances to check our results such as that of a ring, a rectangular loop, a solenoid, and a spiral. We used Mathematica 7 and 9 on a Dell Optiplex 2.33 GHz computer to perform the calculations. We also compare the results to those of FastHenry (FastHenry2 Porting Version 3.32), to understand the advantages and disadvantages of using this method to those of a state-of-the-art type of simulation tool. We also test our results against a simpler version of air core reactors used in the field, a two-layer, two-coiled spiral, which was built to test this method, since a known solution for this configuration does not exist.

Inductance calculations for a ring
A conductor in the shape of a solid ring is illustrated in Figure 1 These coordinates are illustrated in Figure 1. Note that a is the radius of the ring to the center of the conductor, and b is the radius of the conductor. For a uniform current, the term f 1 : f 0 1 can be represented as, and for surface current, the term f 1 : f 0 1 can be represented as, where d represents the Dirac delta function. The variable r is the distance between two volume elements in the conductor, and can be defined as, which, using the formulas in Eqn.
The volume element dV can be represented as,
Because of the high degree of symmetry of the geometry of the ring, Eqn. (1) can be integrated analytically to obtain an exact solution for its inductance and we can check the numerical solution for different ring sizes against this exact solution given by, for a uniform current, and by for surface current [7]. The appropriateness of using the uniform current model or the surface current model can be determined by computing the skin depth d~1= ffiffiffiffiffiffiffiffiffiffiffiffi ffi pf m o s p , where f is the frequency of the current and s is the conductivity. Table 1 and Table 2 describe the results. If neither a surface current or a uniform current appropriately describes the current distribution, a more precise formula for the current densities can be used in Eqns. (1) or (2); however, for this configuration, and for the others described in this study, the percent differences in the results for surface and uniform current distributions are small, less than 10%, so the error in using the uniform or surface distributions for an in-between distribution would be small. The numerical integration was performed using the global adaptive method, the default integration strategy of Mathematica. The match between the exact solutions of Eqns. (12) and (13) and the results produced by this method, shown in Tables 1 and 2, is excellent with the largest error (% deviation from exact solution) being.76% but with most errors less than.01%. The small error that remained is due to the errors near the singularity, where convergence can be difficult. In fact, for almost all the inductance calculations described in this paper, Mathematica produced a warning that the global error in the numerical integration had not decreased monotonically; yet, as  will be shown in this paper, all the results were highly accurate. Convergence could be obtained by lowering the precision of the result (using the PrecisionGoal option) but it was found that the results were usually not as accurate. Timing of the numerical integration was on the order of seconds for the surface current model, and tens of seconds for the uniform current model. FastHenry, which uses rectangular segments, defined by the user, produced slightly less accurate results, with the ring discretized into a 100 segments, with cross-sectional areas set to be the same as those of the rings. However, it only took seconds for both uniform currents (low frequency input) and surface currents (high frequency input), although time was needed to prepare its input files.

Inductance calculations for a rectangular loop with a rectangular cross-section
A conductor in the shape of a rectangular loop (sides of lengths l 1 and l 2 ) with a rectangular cross-section (height, h, and width, w) is shown in Figure 3. From symmetry, we need only compute the self-inductances of two sides with different lengths, and the two mutual inductances between opposite sides and multiply their sum by two to obtain the self-inductance of the loop. The two perpendicular sides do not add mutual inductance because the current densities are perpendicular to each other and therefore f 1 : f 0 1~0 . For the self-inductance of a side, for a uniformlydistributed current, f 1 : f 0 1~1 =(hw) 2 , and for the mutual inductance of the opposite side, f 1 : f 0 1~{ 1=(hw) 2 . In this case, of course, we integrate in rectangular coordinates. To make the integration for the self-inductance more tractable, the numerical integration is performed with the singularities at the limits of the integration range. This geometry also has much symmetry and Eqns. (1) and (2) had been integrated analytically to obtain an exact solution [10], Table 2. Comparison of results between numerical method and exact solution for inductances for rings of surface currents (Eqn. (13)).     Table 3 describes the results. Again an excellent match is achieved between our numerical solution and the exact solution with the largest error being.65%. FastHenry performed slightly more accurately with its largest error being.17%. However, FastHenry was much faster, taking fractions of a second to run whereas this method took tens of seconds, which makes sense for a rectangular configuration, since FastHenry uses rectangular elements, and only four elements need to be specified as input.

Inductance calculations for a solenoid
A conductor in the shape of a solenoid with N turns is illustrated in Figure 4. This calculation is nearly the same as that of the inductance of a ring except that in Eqn. (3), the z coordinate depends on the turn number i, z i~2 bizr sin w, ð15Þ Table 3. Comparison of results between numerical method and exact solution for inductances for a rectangular loop with a rectangular cross-section (Eqn. (14)).   where i varies from 0 to N{1, and the integration is performed over each turn. A first approximation to a solenoid is a cylinder with a surface current, whose inductance, because of its symmetry, can be obtained by analytically integrating Eqn. (1) to obtain Lorenz's exact solution L base [11]. To obtain a better approximation to the solenoid inductance, L sol we subtract Rosa's empirically derived correction factors to account for the effect of the coils [12].
L base is given by, where a~radius of the solenoid in meters; N~number of turns; l~length of the solenoid in meters; w~a=l; and m~4w 2 : K(k) and E(k) are complete elliptic integrals with K(k)~Ð L corr can be written as, where A and B are empirically derived constants; A depends on the ratio of the diameter of the conductor to the turn pitch, and B depends on N. Values for A and B have been tabulated [13] and numerical values for all these quantities can be found online [14]. Table 4 describes the results for solenoids of various radii and number of turns. The match is again excellent with errors less than 1.4%. This method worked considerably better than FastHenry. For 5-turn solenoids, with 100 elements, its largest error was 5.1%. Increasing the number of discretization segments did not reduce the error. This method ran for 1-3 minutes and FastHenry took 2 minutes. For a 100-turn solenoids with 1000 elements, FastHenry ran out of memory and did not produce a result. While this method produced an accurate result for 100-turn solenoids, it did take 1-6 hours to run, with the smaller radius, giving the larger run-time.  Table 5. Comparison of the results of between numerical calculation and Wheeler's approximate formula L Wheeler (Eqn.(20)) for a spiral ( Figure 5).

Inductance calculations for a spiral
Another check on this method is to compare L self for a flat spiral ( Figure 5) with that of Wheeler's empirical formula [15]. To determine L self , we use the coordinates described in Eqn. 3 with a, the radius of the spiraled conductor, described by, where, a o is the widest radius of the spiral (not including the outermost conductor), p is the pitch of the turns which equals the diameter 2b of the conductor plus the spacing between the conductors. Note, for the spiral, there is no axisymmetry and we can no longer let h{whzh' (although we continue to let w{wwzw'); moreover, the limits of h are from 0 to 2pN. The empirical formula is given by, where L Wheeler is the inductance in mH, r is the radius to center of the turns in inches, and w is the width of the turns, in inches. We can obtain r and w using the following formulas, where p is the pitch in inches and D i is the inner diameter of the spiral, which can be obtained from, where the factor 39.37 converts from meters to inches. Wheeler's formula is accurate to within 5% when there are sufficient number of turns (ww:2r), when the spacing between the turns is small, and when the skin effect is not important. As can be seen from Table 5, the results for L self match to within 10% the results of L Wheeler for the uniform current model, and match to within 6% for the surface current model. The match for FastHenry was considerably worse, with the largest error of 42% for a uniform current and with the largest error of 39% for surface current. Doubling the number of elements in the input file for FastHenry did not reduce the error. The run-time for this method was a few minutes and did not increase with the number of spiral turns while for FastHenry, the run-time for the 4-turn spiral, with 100 elements was 20 seconds, and the run-time for the 15-turn spiral which had 500 elements was 10 minutes.
Inductance calculations for a two-layer spiral with two coils Shown in Figure 6 is a two-layer spiral with two coils that was built to test this method. As can be seen from the figure, this configuration consists of two interlaced coils that are arranged in two layers of spirals. From basic circuit theory of non-ideal transformers [16] when a transformer's inductance is measured between the ends of each coil with the coils connected in parallel, its inductance is equal to, By using Eqn. (3) as applied to a single layer spiral (Eqn. (19)), with Eqn. (15) (since we have two layers) for each of the two coils, and also by computing their mutual inductance (Eqn. (2)), we can compute L 2layer . The dimensions and the inductances of both the model and the measurements are shown in Table 6. The Figure 6. A two-layered spiral transformer with two coils, constructed to test the method presented in this paper. Left: The top view; the two spirals are interlaced and there are two layers of the them. Right: Side view. The formula for its shape is described in Eqn. (3), (15), and (19); its dimensions are given in the caption of Table 6. doi:10.1371/journal.pone.0111643.g006 Table 6. Measurements of the two layered spiral shown in Figure 6. The dimensions are a o~: 086m, b~:001m, pitch~:008m, and N~5, and the measured inductance is L measured~4 :4mH. measurements were made using an Agilent LCR meter. We see that the results match the measurements to within 2%, which is likely within the accuracy of the measurement considering stray inductance at the ends of the spiral and of the probe leads.

Discussion
From the results for inductors for which the exact solutions for their inductances are known, such as that of a ring, we see that numerically integrating the inductance formula derived from energy formulations is highly accurate, and likely more accurate than previously used empirical methods that do not claim to give a highly accurate estimate. Moreover, this method, in theory, can be used for any configuration of inductors. The disadvantages of this method are that the length of time of the numerical integration may be long, if the program is run on a slow computer, and errors may be introduced because of the difficulty of convergence near the singularity. We performed the numerical integration inside of Mathematica which meant that it ran on a virtual machine, increasing its run-time, instead of directly on the CPU. Running the code in C directly on the CPU would likely eliminate any disadvantage of long run-times.

Limitations of the study, open questions, and future work
This study showed that inductance calculations based on first principles can be performed on several different coil formations using a sophisticated mathematical computing language which performed the numerical integration of the singularity so that the researcher does not need to become expert in numerical integration techniques or expert in how a simulator such as FastHenry or Ansys Maxwell works. However, the researcher still needs be cognizant of the underlying principles of the inductance calculation which may be important in sharpening his skills in inductor design. Note that this study was limited to only five types of inductor configurations; it stands to reason that this method can be applied to other configurations although this needs to be proved. Another caveat is that almost always the global error of the numerical integration did not decreased monotonically, although accurate results were obtained nonetheless. An additional drawback of this method is that the run-time becomes large when simulating more complex configurations. However, FastHenry did not produce a result for these complex configurations although it is suspected that a greater understanding of its workings may allow the user to fix the problem. Future work might entail performing the numerical integration outside of the Mathematica environment so that the calculation runs directly on the CPU and therefore the runtime would be decreased.

Conclusion
A versatile method based on first principals of computing the inductances of conductors has been described and shown to be accurate, with its accuracy limited by the technique used to numerically integrate the inductance formula. This method allows for the calculation of an inductance using the precise geometry of the inductor particularly an air reactor, or transformer, whose inductance may not have been previously derived or present in look-up tables and empirical formulas. It also eliminates the need to for the user to focus on the numerical aspects of the problem at the expense of the physics of the calculation. Since, in theory, this method will work for any geometry, it can facilitate the design of inductors and transformers to meet particular specifications. This method could not have been used before the advent of fast computers due to the computationally intense method for numerical integration; however, now it is practical and likely will be used to develop inductors and transformers whose inductances would previously haven been exceeding difficult to derive. This method proved to be more robust for calculating the inductances of certain solenoids and spirals for at least one other open-sourced simulator.