Numerical Methods for
Inductance Calculation
Part 2 – Inductance Calculation Refinements
D) The Helical Formula, and Coils of Large Pitch
Introduction
When I began writing Part 2c, my goal was simply to code Snow's inductance formula as a simple computer program function. Snow's formula's main attraction was that it accounted for the helicity of the inductor, while the previous examples discussed on these pages did not. Upon reviewing Snow's derivation however, it became clear that his formula contained certain approximations making it unsuitable for coils of large pitch. This is rather disappointing, because it had seemed initially to be a promising method for analysing such coils. However it was also clear, that the formula could be modified, removing the problematic approximations, making it suitable for coils of any pitch. To do this, it was necessary to re-derive Snow's formula for the mutual inductance of two helical filaments, replacing the independent variable θ (angle of rotation of the filament winding) with a new variable s (distance along the filament). It was also necessary to replace the variable px (coil pitch factor), which could become infinite in certain situations, with the variable ψ (coil pitch angle) which is always finite. The details of the revised formula are discussed below.
Limiting Cases
It's always a good idea to cross check results of different formulae and calculation methods whenever possible. The Snow formula is the only one encountered so far that takes helicity into account. This suggests two limiting cases for verifying the formula. The first is when the pitch angle is 0: a circular conductor. The second is when the pitch angle is 90°: a straight conductor.
For the case of ψ=0, the coil can only be a single turn or a part thereof, (because with no pitch, the second and successive turns of wire would attempt to occupy the same space as the first turn). The case of a single circular conductor was treated in Part 1b, where it was demonstrated that its inductance could be calculated using a formula by Maxwell. Maxwell's formula is in good agreement with Snow's formula for the ψ=0 case.
Now, moving beyond the single turn case, consider the case where we use a fixed length ℓW of wire, say one meter, and a coil form of fixed radius r. We are free to wind a coil with any pitch, winding as many turns as necessary until the wire is used up. In this case then, the number of turns will vary with the pitch angle ψ:
Not surprisingly, as the pitch increases, the number of turns decreases (as we are not stretching the wire). However, as the pitch angle approaches 90°, the "coil" starts to approximate a straight conductor, and when the pitch angle is exactly 90°, then the inductor is indeed a straight conductor whose inductance formula is given by Grover [2], page 35 as:
where ℓW is the wire length and dW is wire diameter. (In this formula, dimensions are in cm, and inductance is µH.)
Unfortunately, the calculation method, as implemented in Part 2c, fails in situations where pitch angle approaches 90° (and is highly suspect for moderate values of ψ), giving an inductance value several orders of magnitude too low. Part of the reason is the simplified term in Snow's formula, as previously discussed, and which is easily corrected. However, there are three other problems which must also be addressed.
Problem 1
Snow's helical mutual inductance formula assumes that the two filaments have the same starting angle of 0, and the same ending angle of 2πN. This is fine for coils of small pitch because the distance between filaments is approximately equal to the displacement of the second filament along the y-axis. However, as the pitch becomes significant, it is obvious that the distance between filaments is considerably less than the offset along the y axis. This is illustrated in the following diagram looking in the x direction at the yz plane.
For a pitch angle ψ=45°, the distance between the filaments g is only cos(45°) or 0.7071 times the offset y0 along the y-axis, and in general, the distance between the filaments is cos(ψ) times the y offset.
It is clear that in order to space the filaments apart by a specified amount, and to maintain a reasonable spatial relationship, the proper solution is to adjust the starting angle of the second filament as well as giving it a y offset as shown below. The combination of starting angle and y offset for the second filament are chosen so that a line connecting the starting points of the two filaments will be perpendicular to the filaments' tangents. Therefore, the mutual inductance formula must be revised to allow filament 2 to have a non-zero starting angle. This angle will be denoted as θ0.
The separation between filaments is denoted in the diagram by the variable g, and has a value given by:
g=√(y02+(rθ0)2)
In practice, the value g will be set equal to the gmd of the conductor, and the values of y0 and θ0 will be calculated so that the line joining the starting points of filament 1 and filament 2 is perpendicular to the tangents to the filaments. In the extreme case when ψ=90°, then y0 will be zero, and filament 2 will simply be rotated in the xz plane.
Problem 2
When the pitch angle ψ approaches 90°, the pitch factor px approaches infinity, and so to handle the case of ψ=90°, the pitch factor cannot directly appear in the formula. Therefore, the formula will be rearranged to use ψ and functions thereof which will always be finite.
Problem 3
Also related to the case when the pitch angle ψ approaches 90°, the angle θ (angle of rotation of the winding), which is used as one of the variables of integration, will approach a constant value. Evaluating an integral between two limits, which are the same value, is not a particularly rewarding experience. Therefore, the formula will be revised by replacing θ with s which is defined here as the distance along the filament, measured from the starting end.
The New Mutual Inductance Formula
The new formula is derived starting with the Neumann Integral, and is somewhat involved. Rather than present all of the gory details on this page, I have given a detailed derivation in reference [21], and will simply give the resulting mutual inductance formula here. The new formula for the mutual inductance of two helical filaments is then:
Where the start of filament 1 is positioned at x=r, y=0, z=0; s1 and s2 refer to distances along filaments 1 and 2 respectively, and where:
ℓW = Total length of the filament
r = radius of winding
kψ = cos(ψ)/r
y0 = Offset of filament 2 along the y-axis
θ0 = Angular offset of filament 2 in xz plane
θ0 and y0 are functions of r, ψ, and the required distance g separating the filaments. The calculation of these values will be discussed later.
Getting formula (23) into the form of a single integral will parallel the method used by Snow, and as described briefly in Part 2c. The variable φ is substituted for (s1-s2). For clarity, the integrand of formula (23) will be represented by the function f(φ) such that:
Then the integration is carried out with respect to s2 first, giving:
This is now in a form that can be evaluated using Simpson's Rule. As before, to convert this from a mutual inductance formula to a self inductance formula, the value g, the distance separating the filaments will be set to the gmd of the conductor. The only remaining obstacle, is to calculate the values of θ0 and y0.
Calculation of θ0 and y0
As mentioned above, θ0 and y0 are functions of r, g, and ψ.
When ψ=0, then y0=g and θ0 =0.
When ψ=90°, then y0=0; and θ0 is approximately equal to g/r (assuming dW <<r). The exact value for θ0 is slightly different due to the curvature of the cylinder on which the coil is wound. The exact value is:
θ0 =arccos(1−g2/(2r2))
When ψ is a value other than 0 or 90°, there is no closed form formula for θ0. But the following expression is exact and can be evaluated numerically:
g2−2r2 +r2(2cos θ0 − θ02 tan2 ψ) = 0
In the program code which follows, this expression will be evaluated using the Newton-Raphson Method. which solves for θ0 in just a couple of iterations when using a starting value θ0 = g/r sin ψ.
Once θ0 has been calculated, then the value for y0 is given by:
y0 = √(g2−2r2(1− cos θ0))
Program Code
The Open Office Basic code for the helical inductance formula is given here:
Function HeliCoilS (byVal Lw as double, psi as double, r as double, dw as double, MaxErr as double) as double
' Version 1.0, 2011-03-25
' Uses helical filament mutual inductance formula
' evaluated using Simpson's rule, and conductor gmd
' Lw = total length of wire
' psi = pitch angle of winding
' r = radius of winding
' dw = wire diameter
' MaxErr = max allowable error
'
' If Lw>2*pi*r, check that pitch angle >= psi-c (close wound pitch)
If Lw>2*pi()*r then
sinpsic=dw/(2*pi()*r)
psic=atn(sinpsic/sqr(1-sinpsic*sinpsic)
if psi<psic then
' pitch angle is too small,
' so set value of function to an illegal value and exit
HeliCoilS=1e200
Exit Function
end if
End if
' gmd of solid round conductor. Other values may be substituted
' for different conductor geometries
g=exp(-.25)*dw/2
rr=r*r
psio=0.5*pi()-psi
' Calculate Filament 2 offset angle
' Trap for psi=0 condition in which case ThetaO=0 and Y0=g
' Trap for psio=0 condition in which case use simplified formula for ThetaO and Y0=0
' which happens with circular (non-helical) filament
if psi=0 then
ThetaO=0
Y0=g
ElseIf psio=0 then
cosThetaO=1-(g*g/(2*rr))
ThetaO=-abs(atn(sqr(1-cosThetaO*cosThetaO)/cosThetaO)
Y0=0
Else
' Use Newton-Raphson method
k1=(g*g)/(2*r*r)-1
k2=tan(psio)
k2=0.5*k2*k2
t1=g/r*sin(psi)
do
t0=t1
t1=t0-(k1+cos(t0)-k2*t0*t0)/(-sin(t0)-2*k2*t0)
Loop until abs(t1-t0)<1e-12
ThetaO=-abs(t1)
' Calculate Filament 2 Y-offset, using formula (29)
Y0=sqr(g*g-2*rr*(1-cos(ThetaO)))
End if
' Psi constants
c2s=cos(psi)^2
ss=sin(psi)
k=cos(psi)/r
' Start of Simpson's rule code
a=0
b=Lw/32768
If b>Lw then b=Lw
grandtotal=0
Do while a<Lw
dx=b-a
m=1
CurrentErr=2*MaxErr
kat=k*a
kbt=k*b
Sum2=(Lw-a)*(Ingrnd(-a,-kat-ThetaO,ss,c2s,rr,Y0)+Ingrnd(a,kat-ThetaO,ss,c2s,rr,Y0)) + (Lw-b)*(Ingrnd(-b,-kbt-ThetaO,ss,c2s,rr,Y0)+Ingrnd(b,kbt-ThetaO,ss,c2s,rr,Y0))
' Initialize LastResult to trapezoidal area for test purposes
LastIntg=Sum2/2*dx
Do While CurrentErr>MaxErr or m<512
m=2*m
dx=dx/2
Sum=0
SumA=0
for i=1 to m step 2
phi=i*dx+a
kpt=k*phi
Sum=Sum+ (Lw-phi)*(Ingrnd(-phi,-kpt-ThetaO,ss,c2s,rr,Y0)+Ingrnd(phi,kpt-ThetaO,ss,c2s,rr,Y0))
next
Integral=(4*(Sum)+Sum2)*dx/3
CurrentErr=Abs((Integral)/(LastIntg)-1)
LastIntg=Integral
Sum2=Sum2+Sum*2
Loop
grandtotal=grandtotal+Integral
a=b
b=b*2
If b>Lw then b=Lw
Loop
HeliCoilS=1e-7*grandtotal
end function
Function Ingrnd (byval phi as double, kphitheta as double, sinpsi as double, cos2psi as double, rr as double, y as double) as double
'Integrand function called by HeliCoilS()
Ingrnd =(1+cos2psi*(cos(kphitheta)-1))/sqr(2*rr*(1-cos(kphitheta)) +(sinpsi*phi-y)^2)
End Function
The implementation is similar to that used in Part 2c. Again, the value of the integrand has a sharp spike near φ=0, making it necessary once again to use unequal intervals, the smallest near φ=0, and doubling in size as φ increases. As before, each of these intervals is subdivided into smaller equal intervals as many times as necessary, in order to get the specified accuracy.
Unlike the previous version in Part 2c, which had inputs of N (number of turns) and p (pitch), these parameters have been replaced with ℓW (total wire length) and ψ (pitch angle), which are not normally specified when winding a coil. This is only a minor inconvenience, because for coils of finite pitch, the values of ℓW and ψ can be calculated from N and p. (I should again point out the distinction between pitch p and pitch factor px, which differ by a factor of 2π.) The advantage, as mentioned above, using parameters ℓW and ψ, allows us to calculate inductance in the limiting case when the wire is pulled into a straight line. For the case where pitch is finite, the input parameters can be calculated from the usual variables as follows:
px = p/(2π) = ℓCOIL/(2πN)
ψ = Arctan(px/(πD)) = Arctan(p/(2π2D))
ℓW = πND/cos ψ
where:
D is coil diameter = 2r;
p is centre to centre turns spacing;
ℓCOIL is length of coil.
One other item not previously discussed, is that the program checks to see that if ℓW is greater than πD then the minimum pitch angle cannot be less than the close-wound pitch ψC, as the winding would overlap. The close-wound pitch angle ψC is given by:
ψC = arcsin(dw/(πD))
If the specified pitch angle ψ is less than ψC, then the calculated inductance is set to an error value, and the program exits.
On the subject of interval sizes, the hard-coded value of 32768 is used in the program line:
b=Lw/32768
to set the first interval to 1/32768th of the total conductor length. Each successive interval is double the size of the previous one. However, each of these intervals will be split in half a fixed minimum number of times, and evaluated in the inner loop, to ensure that the MaxErr condition is satisfied. The value 1/32768 was set quite small to be conservative, and some experimentation may be in order to find the optimum value. However the optimum is known to vary depending on the number of turns of the coil, the more turns, the smaller the increment required. Likewise, the hard-coded minimum iteration count of 512 in the inner While loop may be subject to some optimization.
Results
The point of all this was to see what happens when we take a fixed length of wire, and a coil form of a given diameter which is at least as long as the wire, so that a coil may be wound of any pitch (varying from 0 to ∞ as ψ varies from 0 to 90°), and we can see how the inductance varies with pitch. For an example, a conductor of length 5 meters, and diameter of 1 mm will be used, and coil diameter is fixed at 25 mm. The following graph shows the results:
The blue line is the inductance. As the pitch angle increases, the inductance drops off and then rises again, reaching the straight conductor inductance value (red dashed line), as it should, when ψ=90°. The violet line indicates how the number of turns of the coil decreases as ψ increases, reaching exactly zero when ψ=90°.
The following graph gives a comparison between the current formula, the unmodified Snow formula, and Rosa's formula:
Rosa's formula is in agreement until ψ exceeds about 5°. Snow's original formula agrees up to about ψ =15°.
The following graph shows the same comparison expanded to show only the range ψ=0..12°:
The curves for Snow's formula and the modified helical formula are indistinguishable, while the curve for Rosa's formula is very close. It should be pointed out that a pitch angle of 5° amounts to quite significant pitch, not likely to be encountered in most coils except perhaps in VHF work. Hence, Rosa's formula appears to be completely suitable for most practical coils, and at a great savings of computation.
As the above curves are too cramped to read accurate numerical values, the following actual program output values are provided for comparison:
Coil 1 (as shown in above graphs):
ℓW = 5000 mm; dw = 1 mm; D = 25 mm; ψ = ψC = 0.72953° (close-wound)
LHEL = 32.8527 µH (new helical formula)
LSNOW = 32.8474 µH (Snow's formula)
LROSA = 32.5840 µH (Rosa's formula)
ℓW = 5000 mm; dw = 1 mm; D = 25 mm; ψ = 5°
LHEL = 6.8808 µH
LSNOW = 6.8304 µH
LROSA = 6.6018 µH
ℓW = 5000 mm; dw = 1 mm; D = 25 mm; ψ = 15°
LHEL = 5.0111 µH
LSNOW = 4.6914 µH
LROSA = 3.8823 µH
ℓW = 5000 mm; dw = 1 mm; D = 25 mm; ψ = 90° (straight conductor)
LHEL = 9.1536 µH
LLINEAR = 9.1535 µH (Grover's straight conductor formula)
LSNOW = N/A
LROSA = N/A
Coil 2 (a single turn when ψ =0):
ℓW = 3141.59 mm; dw = 1 mm; D = 1000 mm; ψ = 0
LHEL = 4.5473 µH
LSNOW = 4.5473 µH
LROSA = 4.5473 µH (coil length set equal to dw)
LMAXWELL = 4.5473 µH (Maxwell's loop inductance formula)
ℓW = 3141.59 mm; dw = 1 mm; D = 1000 mm; ψ = 90°
LHEL = 5.4594 µH
LLINEAR = 5.4593 µH
As can be seen, the new formula shows good agreement with the Snow and Rosa formulae when ψ is small, and good agreement with the straight conductor formulae when ψ is 90°.
Conclusion
I should point out that, while the formula presented here agrees with other formulae in cases where such formulae are available and known to be accurate, i.e., small pitch, and infinite pitch, I'm not aware of any formula which may be used as a cross-check in the intermediate pitch values (ψ=45°, for example), and I have not attempted to verify the formula by experiment. Therefore, until further work is done to verify it, this formula should be used with caution.
The program presented here is much more computationally intensive than Rosa's formula, and consequently, it would be of interest to find a way to modify Rosa's formula to extend its range of accuracy. One such modification, which I've looked at only briefly, is to combine Rosa's formula with the straight conductor formula using weighting factors depending on the pitch angle. Rosa's formula would have a weight of 1 at ψ=0 and a weight of 0 at ψ=90°, while the weighting factor for the straight line formula would be the reverse. The weighting factor could vary linearly with ψ, or sinusoidally, or as a power law. These all show some promise, but more work needs to be done.
Continue to:
Part 3 – Optimizers, Solvers and Empirical Methods
Back to:
Numerical Methods Introduction
This page last updated: February 25, 2021
Copyright 2011, 2015, Robert Weaver