Part 2 – Inductance Calculation Refinements

100. C)Snow's Helical Inductance Calculation Method
     

Introduction

The calculation to be described here is not really a refinement, but a totally different calculation method. However, the discussion fits in well here, because it does follow from the discussion of helical corrections in the previous section. It also introduces another useful computational tool.

Snow's Formula

In 1926, Chester Snow of the US National Bureau of Standards published a formula for the inductance of a helix [19]. Existing inductance formulae, though sufficiently accurate for virtually all practical purposes, ignored the helicity of the winding, and there was no simple way to tell how much effect these neglected components had on coil inductance. With improved technology for measuring inductance, Snow felt there was a need to develop a formula which would fully account for all aspects of the coil geometry. This was an ambitious undertaking involving ninety pages of dense  mathematics.

Snow defined the coil geometry in cylindrical co-ordinates. Rather than representing a point in space in the form (x,y,z), it is represented as (r,θ,y) where r and θ are standard polar coordinates radius and angle, and y is the displacement into the remaining dimension. Consequently, a helix may easily be defined by it's radius r, starting angle θ1 ending angle θ2, and either its pitch or its total length.

The foundation of the work was the development of a formula for the mutual inductance of two helical filaments of N turns, radii r1, r2, axial offset y1, y2: and pitch 2πpx.

where

Note that Snow defines the pitch—the axial spacing of turns—as 2πp in his paper, while p is used for pitch everywhere else on this site. Therefore I've included the x subscript and used px in his formula such that p= 2π px, in order to avoid confusion. The variable px will be referred to as pitch factor rather than pitch.

The mutual inductance formula is a straightforward derivation from the fundamental principles of magnetics and three dimensional geometry. In working through the derivation of formula (20) however, I noted that Snow seems to have simplified a term. It appears that the numerator of the fraction outside of the integrals should be

rather than simply r1r2. Whether this was intentional or not is unknown. Its omission has almost no effect on coils of small pitch, and Snow comments in a follow-up paper that he is primarily concerned with coils of the type of geometry used by the Bureau of Standards for their standard inductors which were typically of small pitch. So, it is possible that he deliberately left the term out. However, he does imply in paper 537 that the mutual inductance formula (20) is completely general, which does not seem to be the case. For the time being we will leave formula (20) as Snow published it, but we will revisit this later as we look at coils with larger pitch.

While the derivation of formula (20) is straightforward, the derivation of an analytical solution to this integral is a very messy process which accounts for about half of the mathematics in the paper.

Snow makes the substitution φ=(θ2-θ1) so that:

Then integrating first with respect to dθ1 gives:

where, now:

So far, so good. But now, things start to get ugly. There is no simple analytical solution to these integrals. Dr. Snow uses a number of transforms and series expansions to arrive at—several dozen pages later—a closed form expression for the mutual inductance. Even so, he applies numerous simplifications, dropping terms greater than second order. Getting a closed form expression for M is important because, in order proceed to get an expression for self-inductance, it is still necessary to integrate this helical filament mutual inductance expression over the cross section of the conductor—a further double integration.

At this point we will part ways with Snow's derivation, and take a different approach. We have, in integral form, formula (21) for the mutual inductance of two helical filaments. In Part 1b we used a filamentary mutual inductance formula along with the geometric mean distance of a conductor cross-section to calculate the self-inductance of a coil (as a series of parallel coaxial rings). We will take the same approach here, using the geometric mean distance of a round (or thin tubular) conductor and the mutual inductance formula to calculate the self-inductance of a helical conductor of finite cross-section. We will use Simpson's rule to evaluate the integral numerically.

Simpson's Rule

Simpson's Rule is a numerical method for evaluating integrals. To evaluate the integral of a function f(x) over the range x=[a..b], the function to be integrated is evaluated at three points a, b and m, where m=(a+b)/2, the midpoint between a and b. The area under a parabolic curve which passes through the three points f(a), f(b) and f(m) is used as an approximation of the true value of the integral.

The approximate value of the integral, using the parabolic curve, is then given by:

It should be obvious that the accuracy of this approximation will depend on how well a parabola fits the actual function. In practice, the interval over which we wish to integrate will be broken into smaller sub-intervals, and Simpson's Rule is applied to each sub-interval. If we break the interval into n subintervals, then Simpson's Rule gives:

A few things to note:

1. •There is always an even number of intervals, which means that there will be an odd number of points evaluated;
   

2. •The odd numbered points have a coefficient of 4;
   

3. •The even numbered points (except x0 and xn) have a coefficient of 2;
   

4. •The end points x0 and xn have a coefficient of 1;
   

5. •The more points evaluated, the more accurate the result.
   

It's apparent that we can evaluate the function for all of the values of x using a program loop, and apply the correct coefficient depending on whether the subscript of the x value is even or odd. However, a more efficient method is to maintain two summations, one for the even values and one for the odd values, and then multiply these sums by the appropriate coefficient. A less obvious reason for doing this is that if using n sub-intervals doesn't give us the desired accuracy, we can double the number of sub-intervals, and the points that had been previously calculated won't have to be discarded. All of the previous odd and even points will become the new even points, and only the new odd numbered points (which will fall between the old points) will be evaluated.

The following Open Office Basic code uses this method to implement Simpson's Rule:

Function Simpson (ByVal a as double,b as double,MaxErr as double) As Double

  'Numerical integration using Simpson's Rule

  'Integrates the function 'TestFn(x)' between limits x=a and x=b

  'Note that TestFn must be defined elsewhere as another function

  'Stops when the difference between subsequent calculations

  'is less than MaxErr

  dx=b-a

  n=1

  CurrentErr=2*MaxErr

  Sum2=TestFn(a)+TestFn(b)

  'Initialize LastResult to trapezoidal area for test purposes

  LastResult=Sum2/2*dx

  Do While CurrentErr>MaxErr

    'Double the number of points n each time through the outer loop

    'Given the numbering convention that point Xo=a (lower limit)

    'and Xn=b (upper limit), then only the odd numbered points need

    'to be calculated. The even points are all of the points

    'calculated from previous iterations (beginning with the two

    'endpoints calculated outside of the loop),

    'and stored in variable Sum2

    n=2*n

    dx=dx/2

    Sum=0

    for i=1 to n step 2

      Sum=Sum+TestFn(i*dx+a)

    next

    CurrentResult=(4*Sum+Sum2)*dx/3

    CurrentErr=Abs(CurrentResult/LastResult-1)

    LastResult=CurrentResult

    Sum2=Sum2+Sum*2

  Loop

  'msgbox str(n+1)+" Function Evaluations"

  Simpson=CurrentResult

End Function

Function TestFn (ByVal x as Double) as double

  'Used for testing Simpson's Rule function

  'The following is just an arbitrary function

  TestFn=.25*x^3+9/8*x^2+3*x

end Function

The Simpson's Rule function is heavily commented for clarity. The actual executable code is quite short. The function to be integrated is defined as a separate function TestFn() so that it may be called from several different places within the main routine. The function starts with a single interval. The end points are calculated outside of the loop, and the single midpoint is calculated inside the for- next loop which does only a single iteration the first time around. Each time through the outer loop the number of sub-intervals is doubled.

The termination test compares the last calculated value of the integral with the current value. If the relative difference is less than the specified value, MaxErr, then the loop terminates.

Simpson's Rule Applied to Inductance Calculation

Returning to the problem at hand, using the above format we can now look at evaluating the two integrals of mutual inductance formula (21). Before jumping in though, there are a few simplifications which can be made.

For this calculation, the two filamental radii will be equal. Hence:

r1=r2=r

We will set y1=0 and y2=g, the geometric mean distance of the conductor cross-section.

We must also multiply the entire formula by µ0/(4π) = 10-7 in order to get a result in the SI unit of Henries (with input parameters in meters). Snow's original formula gives inductance in units of length, a system of measurements which was not uncommon at the time, but is now archaic.

Formula (21) then becomes:

and:

We thus use LH to denote the self-inductance of a helical conductor.

The geometric mean distance (gmd) of the cross-section of the conductor is almost the same as the value used in Part 1b when we were dealing with Maxwell's coaxial circles formula. However, there is a slight difference. The cross-section used in this formula is parallel to the y-axis which is not normal to the direction of conductor. It's the same situation as described in Part 2b where cutting the conductor parallel to the axis of the winding results in an elliptical cross section. So, we need to use a gmd formula for an ellipse rather than a circle. This formula is given in Grover [2], page 21 as:

where a and b are the semi-major and semi-minor axes of the ellipse. The semi-minor axis will simply be the conductor radius, while the semi-major axis will be equal to the conductor radius divided by the cosine of the pitch angle ψ where:

ψ=arctan(p/πd)

(d is coil diameter, p is winding pitch)

Then the gmd is given as:

for a solid conductor at low frequencies (dw is wire conductor diameter).

Or, for a thin tubular conductor, or a solid conductor at high frequencies:

Finally we are at a stage to present a routine to calculate the self-inductance of a helical conductor. Written in Open Office Basic, it is:

Function HelicalCoil_1 (byVal N as double, pitch as double, r as double, dw as double, MaxErr as double) as double

'Uses Snow's helical filament mutual inductance formula

'evaluated using Simpson's rule, and conductor gmd

  p=pitch/2/pi()

  psi=atn(pitch/(2*r*pi())

  g=.25*(1+1/cos(psi))*exp(-.25)*dw 'gmd of conductor

  Ntp=2*pi()*N

  rr=r*r

  a=0

  b=Ntp

    dx=b-a

    m=1

    CurrentErr=2*MaxErr

    Sum2=HlxIntgrnd(1,-a,Ntp,p,rr,g)+HlxIntgrnd(1,-b,Ntp,p,rr,g)+HlxIntgrnd(-1,a,Ntp,p,rr,g)+HlxIntgrnd(-1,b,Ntp,p,rr,g)

    'Initialize LastResult to trapezoidal area for test purposes

    LastIntg=Sum2/2*dx

    Do While CurrentErr>MaxErr

      m=2*m

      dx=dx/2

      Sum=0

      for i=1 to m step 2

        Sum=Sum+HlxIntgrnd(1,-i*dx-a,Ntp,p,rr,g)+HlxIntgrnd(-1,i*dx+a,Ntp,p,rr,g)

      next

      Integral=(4*Sum+Sum2)*dx/3

      CurrentErr=Abs((Integral)/(LastIntg)-1)

      LastIntg=Integral

      Sum2=Sum2+Sum*2

    Loop

    grandtotal=grandtotal+Integral

  HelicalCoil_1=1e-7*rr/(1+p*p/rr)* Integral

end function

Function HlxIntgrnd (byval sn, phi as double, N as double, p as double, rr as double, y as double) as double

'Integrand function called by HelicalCoil()

  HlxIntgrnd=(N+sn*phi)*(cos(phi)+p*p/rr)/sqr(  (2*rr*(1-cos(phi))) +(y+p*phi)^2)

End Function

Now for a few comments. It can be seen that even though we have two integrals to evaluate, there is only a single routine. Both integrals are evaluated at the same time in the same program loop and summed together. It is convenient to do this since both integrals are evaluated over the same interval, the only difference being the sign of one of the limits. The order of calculation of the second integral has been reversed as it makes no difference to the calculation result, but simplifies the program. Only a single function is used for the integrand, as the only difference is in the sign of φ in the first term in the numerator. This is handled by passing a sign parameter in the form of the variable sn which is set equal to either +1 or −1.

As can be seen, there is a very sharp spike near zero, but this spike accounts for a significant portion of the value of the integral. In order to evaluate it accurately, the spike needs to be divided into a sufficient number of sub-intervals. But, in the process of doing this, the remainder of the curve is also unnecessarily divided into a huge number of sub-intervals. A possible solution is to add an outer loop to the program, to break the integration up into unequal sections. The first section would be the interval from 0 to π/2, then from π/2 to π, then from π to 2π, doubling the interval each time through the outer loop. The inner loop would then equally subdivide these sections as many times as necessary to achieve the desired accuracy. Another problem is that because the two integrals are summed together as they are calculated, it creates a problem with the error calculation and termination test. I had originally assumed that since the curves of the two integrands are nearly identical, combining the integrals would not affect the error test. However, on closer examination, there are many cases where the two curves diverge significantly near φ=0. Therefore, I separated the integral summations.

The modified routine follows:

Function HelicalCoil (byVal N as double, pitch as double, r as double, dw as double, MaxErr as double) as double

'Uses Snow's helical filament mutual inductance formula evaluated using Simpson's rule, and conductor gmd

  p=pitch/2/pi()

  psi=atn(pitch/(2*r*pi())

  g=.25*(1+1/cos(psi))*exp(-.25)*dw 'gmd of conductor

  Ntp=2*pi()*N

  rr=r*r

  a=0

  b=pi()/2

  If b>Ntp then b=Ntp

  grandtotal=0

  mtot=0

  'Outer loop

  Do while a<Ntp

    dx=b-a

    m=1

    CurrentErr=2*MaxErr

    SumA2=HlxIntgrnd(1,-a,Ntp,p,rr,g)+HlxIntgrnd(1,-b,Ntp,p,rr,g)

    SumB2=HlxIntgrnd(-1,a,Ntp,p,rr,g)+HlxIntgrnd(-1,b,Ntp,p,rr,g)

    'Initialize Last Result to trapezoidal area for first termination test

    LastIntgA=SumA2/2*dx

    LastIntgB=SumB2/2*dx

    Do While CurrentErr>MaxErr

      m=2*m

      dx=dx/2

      SumA=0

      SumB=0

      for i=1 to m step 2

        SumA=SumA+HlxIntgrnd(1,-i*dx-a,Ntp,p,rr,g)

        SumB=SumB+HlxIntgrnd(-1,i*dx+a,Ntp,p,rr,g)

      next

      IntegralA=(4*SumA+SumA2)*dx/3

      IntegralB=(4*SumB+SumB2)*dx/3

      CurrentErr=Abs((IntegralA)/(LastIntgA)-1)+Abs((IntegralB)/(LastIntgB)-1)

      LastIntgA=IntegralA

      LastIntgB=IntegralB

      SumA2=SumA2+SumA*2

      SumB2=SumB2+SumB*2

    Loop

    grandtotal=grandtotal+IntegralA+IntegralB

    a=b

    b=b*2

    If b>Ntp then b=Ntp

    mtot=mtot+m+1

  Loop

  HelicalCoil=1e-7*rr/(1+p*p/rr)*grandtotal

end function

Function HlxIntgrnd (byval sn, phi as double, N as double, p as double, rr as double, y as double) as double

'Integrand function called by HelicalCoil()

  HlxIntgrnd=(N+sn*phi)*(cos(phi)+p*p/rr)/sqr(  (2*rr*(1-cos(phi))) +(y+p*phi)^2)

End Function

Just adding the few extra steps to incorporate the outer program loop makes a dramatic difference in execution time. The number of function evaluations drops in some cases from 250,000 to 500.

Now that the routine runs at an acceptable speed, we now need to comment on the accuracy of this method of inductance calculation. In initial testing, the values produced by this routine agree with the Lorenz/Rosa method and the Maxwell summation method to within fractions of one percent. This routine tends to produce values slightly higher than the other two methods, especially in cases where the pitch is exaggerated. This is to be expected since the other calculation methods do not account for the axial component of current flow which would contribute to the total inductance.

It's also worth comparing computational efficiency between the various calculation methods so far presented. The Rosa/Lorenz method remains the most efficient, generally requiring about seven or eight iterations of simple arithmetic functions to solve the elliptic integral portion, followed by a couple of single line closed form expressions for the round wire corrections. The Maxwell summation method solves the elliptical integral once for each turn in the coil. So, for a coil of N turns it would iterate approximately 8N times. For coils of a modest number of turns, say 100 or fewer, it's reasonably fast, and should be faster than the Helical Simpson's Rule method. The Helical method requires the most iterations, and each iteration involves two cosines and a square root. From initial testing, the required number of iterations is dependent on the number of turns in the coil. For ten significant figures, it requires from 30 to 50 iterations per turn when the number of turns is less than 100.  Of course, coils with significant pitch are unlikely to have a large number of turns, making this formula acceptably fast in situations where it would be most useful.

Comparison of Formulae

Following is a comparison of the formulae for the same twelve example coils considered in Part 2a of this discussion. The Snow numerical integral calculation is shown as L-Snow. In addition a variation of the Maxwell Summation method is also included. It is designated L-SummationC, and includes the turn length adjustment discussed in Part 2b and a turn spacing correction (also briefly discussed in Part 2b). Relative Error is relative to the L-Lorenz/Rosa value.

In this first set of coils 1 through 4, where the pitch is moderate, the Snow values are noticeably higher than the others, It is interesting to note that the relative difference between Lorenz/Rosa and Snow is a fairly consistent 3500 parts/million. This higher value for the Snow formula is to be expected due to the contribution of the axial current component which is not included in the other formulae. Even so, a difference of about 3500 parts/million is not huge.

In coils 5 through 9, the pitch is greater, and now we see that the relative difference is no longer consistent. It increases with the number of turns in the coil.

In coils 9 through 12, with larger coil diameter, there is now a negative correlation between number of turns and the relative error. We can conclude then, there is no simple relationship correlating the values calculated by the Lorenz/Rosa and Snow methods.

With this limited set of example coils, the only thing we can really say is that the values given by the Snow formula are higher than the Lorenz/Rosa formula by a modest amount. We can also see that the Summation formula with the included helical corrections appear to correct the inductance in the right direction (assuming the Snow formula is the correct one), but don't adjust the values far enough. A further refinement of both the Lorenz/Rosa and Summation formulae might be to add in the inductance of a tubular axial current of the same diameter and length of the coil. This will be addressed in the future.

Summary

On this page, a third method of inductance calculation has been presented. The present method accounts for helicity of the windings. The results are in close agreement with the previously presented calculation methods. It is interesting—and heartening—to see how completely different calculation methods produce numbers which are in agreement with each other.

Snow released a simplified version of his formula in a later (1932) paper [20]. In it, he claimed that the simplified formula was just as accurate as the original. This is quite likely true, as he had had nearly six years to review the work, and no doubt make refinements. However, it's wise not to take things at face value. In at least one place in his discussions, he points out simplifications that have been made because they would not have any significant effect on NBS standard inductors since they are generally all in the millihenry range and small pitch. And, on at least one occasion, we have seen the result of one of these simplifications (and we will see more in the next section). Consequently, one needs to review the approximations and simplifications carefully to ensure that they don't affect the calculation accuracy where the simplifications do not apply.

In the next section, the method is extended to account for shortcomings in the formula when applied to coils of large pitch.

Continue to:

Part 2d – The Helical Formula, and Coils of Large Pitch

Back to:

Numerical Methods 2b

Numerical Methods Introduction

Radio Theory

Home

This page last updated: February 25, 2021

Copyright 2010, 2015, Robert Weaver