Phil Lucht Math & Physics Archive
Home / Math and Physics Files / Physics / Transmission Lines / Notes By Chapter and Appendix / Appendix P eddy currents / Obsolete

Appendix P Eddy Currents

DOCX · 384.0 KB
Open DOCX file

Phil's draft appendix (dated 3.26.05, in an Obsolete folder) for his transmission line notes. It introduces eddy currents via Faraday's law, then a perturbation expansion in a small parameter valid when the skin depth is much larger than the object size, and notes eddy current testing. It works out self-induced eddy currents in a round wire; several sections (thin plates, proximity effect) are only listed in the outline or are not yet written.

AI-written summary; may contain errors.

Extracted text (machine-read; may contain errors)
This is the Title PhL 3.26.05 Appendix P: Eddy Currents and the Proximity Effect 1 P.1 Eddy Current Analysis 1 P.2. Eddy currents in a thin round plate in a uniform B field 3 P.3. Eddy currents in a thin round plate in a B field with a gradient 3 P.4 Self-Induced Eddy Currents in a round wire 3 P.5 The Eddy Currents induced in a quiet round-wire by an external B field 7 P.6 Eddy Currents induced in an current-carrying wire by an external B field 9 P.7 Summary of Round Wire Examples 10 P.8 Eddy currents in Transmission Lines: The Proximity Effect 10 P.9 Quantitative Evaluation of Eddy Currents and The Proximity Effect 12 P.10 Proximity Effect and Wire Resistance 13 Appendix P: Eddy Currents and the Proximity Effect The Maxwell curl E equation (1.1.2) and its integral form are, curl E = - ∂tB C E ds = -∂t[∫S B dS] . (1.1.36) The integral form is often written as Eemf = -∂t[magnetic flux] and one says that a changing magnetic flux through a loop induces a voltage Eemf (an "electro motive force") in that loop which then drives a current around the loop if the loop lies in a conducting medium. This is Faraday's Law of Induction and the loop of interest is usually a thin wire or coil of such wires inside, say, an electric generator. The wire or coil of wires is attached to some Rload and some current I flows through the loop and load. There is Ohmic loss I2Rloop in the generating loop(s), but if Rload >> Rloop this loss is minimal in the context of the generator. When the loop lies inside an open conducting medium, things become more complicated and the currents which are then driven around mathematical loops in that medium are called eddy currents. The word eddy suggests the way water swirls around in a constrained environment when driven by wind or water currents (see plots below). Just as the water flow velocity can have no normal component at a boundary (a steep river bank for example), an electrical eddy current density generally has no normal component at a boundary of the conductor. An exception to this rule occurs if the eddy current is feeding a charge density on the outer surface of that boundary, and this exception would apply to water flow as well if a bank were shallow and can act as a temporary reservoir of water. The analogy is not exact, but the word eddy is apt. In practical terms, eddy currents are normally seen as undesirable, as in a transformer core, since they represent Ohmic loss which results in power waste and heating of the core (Rload = 0). Sometimes, however, eddy currents are useful, such as in non-destructive testing for internal cracks in metal parts, as noted below. In this Appendix, we explore the nature of eddy currents and compute these currents for some simple situations. We then show how one can interpret both the skin effect and the proximity effect (non-uniform current densities in nearby conductors) in terms of eddy currents. P.1 Eddy Current Analysis [ I plan to replace this section with material in Conduction Eddy and Lines doc. ] We start with some "external apparatus" in which current density Jext flows through some wires creating magnetic field Bext everywhere in space. Both Jext and Bext are of the ω domain form F(x,ω), indicating ejωt time dependence. We identify a region of space close to, but disjoint from, the region of the external apparatus where we shall be placing a conductive Device Under Test (DUT). For example, the DUT could be an airplane wheel or it could be a transformer core very intimate with Jext but still disjoint from it. Assume that the DUT is not yet present and we just have the Maxwell curl B equation (1.1.1), curl (Bext) = μ(Jext + Jdisp) (1) which is valid in all regions of space including that in which the DUT will be placed. Here Jdisp = ∂tD = ε ∂tE is the displacement current. Since the current Jext doesn't exist in the region where the DUT will be, the above becomes curl (Bext) = μ Jdisp // in the DUT region, DUT not present (2a) We now cause the DUT to materialize in its prescribed region of space. To simplify things, we assume that the DUT and its surroundings all have the same magnetic permeability μ. This means there is no magnetic boundary at the surface of the DUT so that the field Bext goes right through the DUT surface unaffected. If this is not the case, one must compute the adjusted Bext inside the DUT. The above equation now reads, curl (Bext + Beddy) = μ (Jdisp + Jeddy) // in the DUT region, DUT present (2b) where we allow that the eddy current in the DUT creates its own magnetic field Beddy. If the conductor is a good conductor like copper, we shall assume that we can ignore the displacement current inside the DUT as discussed below (2.2.2). We can check this assumption later in our examples. Then, curl(Bext + Beddy) = μ (Jeddy) . // in the DUT region, DUT present (2c) Meanwhile, the curl E Maxwell equation says curl(Jext + Jeddy) = -jωσ(Bext + Beddy) This Bext within the DUT creates an eddy current according to Faraday's Law (using J = σE) curl (Jeddy) = - jωσ(Bext) . (3) We now define B(0) ≡ Bext (4) α J(1) ≡ Jeddy (5) where α is a dimensionless "scale parameter". We assume that can see from (3) that if ω is sufficiently small, Jeddy will be small, and we make this fact explicit by writing Jeddy as [α J(1)] where J(1) can be "normal size" and then the eddy current smallness (if it is small) is accounted for by α being small. We then rewrite (3) in terms of these new symbols, curl (α J(1)) = - jωσ(B(0)) . (6) These eddy currents [α J(1)] in turn create a new magnetic field [α B(1)] according to Ampere's Law, curl (αB(1)) = μ(αJ(1)) . (7) We then have to account for this new B field on the right side of (6). But that new B field will in turn induce a new and very small current on the left of (6) we will call α2 J(2). We then have curl (α J(1) + α2 J(2) ) = - jωσ(B(0) + α B(1)) . (8) We can continue this iterative description to obtain curl (Σn=1∞ αn J(n) ) = - jωσ(Σn=0∞ αn B(n)) . (9) Thus, we have an infinite expansion for both the current J and the field B inside the DUT in terms of a dimensionless perturbation-theory expansion parameter α, where all the field contributions J(n) and B(n) are normalized to some reasonable range of values. These expansions are then only meaningful if α < 1, otherwise they diverge. So basically our perturbation theory approach is valid only if ω is small enough to obtain α < 1. If ω is small so α <<1, then we may keep just the first terms in the two expansions so that curl (α J(1) ) ≈ - jωσ(B(0)) (10) curl (αB(1)) ≈ μ(αJ(1)) (11) Defining Beddy ≡ αB(1) we restore our former notation to express these equations as curl (Jeddy) ≈ - jωσ(Bext) // time domain: curl (Jeddy) ≈ - σext (12) curl (Beddy) ≈ μ(Jeddy) . (13) The point of the above discussion is to show that (12) is a reasonable approximate equation if frequency ω is low enough such that | Beddy | << | Bext |. In a very ballpark sense, if the DUT has a characteristic size of D meters, each curl operation in the above equations creates a factor of 1/D, so (1/D) Jeddy ~ - jωσ(Bext) => Jeddy ~ - D jωσ(Bext) (1/D) Beddy ~ μ(Jeddy) => Beddy ~ D μ(Jeddy) => Beddy ~ D μ [- D jωσ(Bext)] = -D2j μωσ Bext and then our condition can be approximately written D2 μωσ << 1. Since δ ≡ is the skin depth for the DUT, the condition can be written D2 (2/δ2) << 1 or δ2 >> 2D2. Thus, we expect an eddy current analysis based on equation (12) to be justified if ω is low enough that δ2 >> 2D2 for a particular DUT situation. A better estimate requires actually computing Beddy from (13) and comparing the result to Bext. When α is not small, there are still eddy currents, but they cannot be treated in the perturbation method just presented. One is faced in this case with the two Maxwell curl equations curl J = - jωσ B J = Jext + Jeddy curl B = μ J B = Bext + Beddy which can be combined into the vector Helmholtz equation as shown in (1.5.27) to give (2+β2)B = 0 β2 = -jωσμ One must then solve this equation for B (using appropriate boundary conditions) and then one obtains the current density J from curl B = μJ . The difference J - Jext can then be interpreted as the eddy current. In typical Eddy Current Testing (ECT) systems, the frequency used might range from 10Hz to 1500 Hz. The idea of an ECT system is to try to detect Beddy using a sensitive Hall Effect or SQUID device, and take note of the field pattern produced by a DUT which is "known good" (has no internal cracks in the metal). An internal crack in a bad DUT will alter Jeddy in some way, which in turn causes an alteration in Beddy which can hopefully be detected. Due to the skin depth penetration issue, the useful depth of such non-destructive testing systems might be up to 15 mm (ballpark). Higher ω generates a larger signal, gives more accuracy on the defect size and location, but penetration depth is less, so there is always a compromise. Often scans at different ω values are optimal for different depths of the defect. ECT is a subject of much current interest and many papers have been and are being written. P.2. Eddy currents in a thin round plate in a uniform B field P.3. Eddy currents in a thin round plate in a B field with a gradient P.4 Self-Induced Eddy Currents in a round wire As a prototype example, we consider our well-studied axially symmetric radius-a round wire of Chapter 2. Now the "external apparatus" and the "device under test" are one in the same! In the zeroth order of the above described perturbation theory (very low ω), the current density in the wire is Jext(x,ω) which is perfectly uniform across the wire cross section and flows in the direction. This current density creates a magnetic field Bext in the direction which is obtained from Ampere's Law, 2πr Bext(r) = μ Ienc = μ Jext πa2 = πr2μJext => Bext(r) = (1/2) μ r Jext . In this problem the "external" current density Jext is generated by the wire itself, as if it were somehow its own "external apparatus". A better notation would be Jdc since this is the ω = 0 current distribution, but we continue to use Jext to maintain contact with the previous Section. Similarly, Bext = Bdc . If the skin depth δ is large compared to the wire radius a, we expect the eddy current analysis of the previous section to be viable, and we find that curl (Jeddy) ≈ - jωσ(Bext) where Bext = (1/2) μ r Jext = Bext(r) . In cylindrical coordinates one writes for an arbitrary vector field F, curl F = [ r-1∂θFz - ∂zFθ] + [∂zFr - ∂rFz] + [ r-1∂r(rFθ) - r-1∂θFr ] For a vector field F which is a function only of r one finds curl F = [- ∂rFz] + [ r-1∂r(rFθ) ] Thus (*) becomes these two equations, r-1∂r[r(Jeddy)θ] = 0 - ∂r(Jeddy)z = -jωσ Bext(r) = -jωσ (1/2) μ r Jext The first equation may be written as ∂r[r(Jeddy)θ] = 0 or [r(Jeddy)θ] = C1 or (Jeddy)θ(r) = C1/r from which we must conclude that C1 = 0 and then (Jeddy)θ(r) = 0, so there is no azimuthal eddy current in the wire. The second equation may be integrated from r=0 to r=r to obtain (Jeddy)z(r) - (Jeddy)z(0) = jωσ (1/2) μ Jext !Syntax Error, Idr' r' = jωσ (1/4) μ Jext r2 . Consider now a tiny circular math loop centered on the round wire axis. Since B → 0 as r→ 0, this loop has no flux, so the induced eddy current around this loop → 0. Thus (Jeddy)z(0) = 0 and we have (Jeddy)z(r) = jωσμ (1/4) μ Jext r2 = - (β2/4) Jext r2 where we use the complex wavenumber symbol β2 = - jωμσ from (2.2.4). Thus, the total current density in the wire obtained from eddy current analysis is in the z direction and is given by Jz = Jext + (Jeddy)z = Jext [ 1 - (β2/4)r2] = Jext [ 1 + j ωσμ (1/4)r2 ] Notice that the eddy current contribution is π/2 out of phase with Jext. Since we have assumed δ >> a, it follows that |βa| << 1 ( since β2 = -2j/δ2 ) and thus (|β|2/4)r2 << 1 so the eddy current contribution is very small, as required to use the first term in the perturbation expansion. Writing κ ≡ ωσμ (1/4)r2 = |β|2 (1/4) r2 << 1 we have |Jz| = | Jext | ≈ | Jext | ≈ | Jext | [ 1 + (1/2)κ2] so that = 1 + {|β|2(1/4)r2}2 = 1 + {(2/δ2) (1/4)r2}2 = 1 + {(1/δ2) (1/2)r2}2 = 1 + (r/δ)4 which then exhibits a very slight skin effect and has the same r dependence as (2.3.10), see Fig 2.7. In Chapter 2 we found in (2.2.30) the following exact result for Jz in a round wire operating at ω, Jz(r) = β . (2.2.30) For small ω (small β) one has for small arguments [ Spiegel 24.5 and 24.6 ] J0(x) ≈ 1 - x2/4 J1(x) ≈ (x/2)(1 - x2/8) 1/J1(x) ≈ (2/x) (1 + x2/8) ≈ (2/x) so ≈ (2/βa) (1-β2r2/4) and then Jz(r) = [(2/βa) (1-β2r2/4)] β = [(2/a) (1-β2r2/4)] = [ (1-β2r2/4)] = Jext (1-β2r2/4) in agreement with our eddy current analysis result ***. At higher frequencies where we no longer have δ >> a, the eddy current perturbation expansion diverges and becomes meaningless and one must instead solve the Helmholtz equation stated at the end of the previous section. In Chapter 2 this task was in essence carried out and the skin effect was observed. One can then interpret the skin effect by saying that the eddy currents cancel the DC current density in the interior of the round wire, allowing a net current to exist only at the periphery. In other words, the skin effect is caused by eddy currents. But this is just a manner of speaking, and is like saying that the skin effect is "caused by Maxwell's Equations", which it is. For a moderate skin effect, we can illustrate the eddy currents by crudely plotting them just in the gray plane of the following drawing : Theses qualitative-only plots are for some particular instant in time. We know that the phase of J varies as shown in Fig 2.8, so we attempt to illustrate only the real parts the currents: Re {Jext} (a) Re{Jeddy} (b) Re{Jext+Jeddy} (c) The closed red curves in the middle drawing represent the Jeddy field lines, and these then represent the actual induced "eddies" of current. One could write Jeddy = σ Eeddy and then they are electric field lines. The lines close on themselves because they have no sources: inside the wire ρ = 0 so div Jeddy = 0 and div Eeddy = 0. Recall that a field in general does not have a constant magnitude along a field line. In the geometry of a round wire, the field lines in fact loop around at the ends of the wire. P.5 The Eddy Currents induced in a quiet round-wire by an external B field In this example, we start with our Device Under Test (DUT) which is a straight round wire which carries no current. Some "external apparatus" creates a time-changing magnetic field Bext = Bext as shown in the figure below, where Bext(x,ω) is for the moment constant in space. At the instant in time shown, the time-domain field Bext(x,t) is increasing in the - direction so that - ext(x,t) points in the + direction, out of the plane of paper. The time-domain eddy current equation curl Jeddy = σ [- ext] C Jeddy ds = σ ∫S [- ext] dS . (1.1.36) implies that the flux change through any math loop in the gray rectangle is positive at our time instant, so according to the right hand rule, the eddy currents in the gray rectangle have the following general shape, (d) Fig **. Eddy Currents in a quiet wire (I = 0) induced by a uniform external B field. To justify the linear variation with x, if we assume that away from the ends of the wire nothing varies with z, and if we write Jeddy as J we find that curl J = (∂yJz - ∂zJy) + (∂zJx - ∂xJz) + (∂xJy - ∂yJx) = - σ ext or (∂yJz) + ( - ∂xJz) + (∂xJy - ∂yJx) = - σ ext which produces the three equations ∂xJz = σ ext => Jz(x) = σext x // linear in x ∂yJz = 0 => Jz = Jz(x) only ∂xJy - ∂yJx = 0 satisfied if Jx and Jy = 0 The eddy pattern in the round wire would have the same general appearance in any slice of the wire parallel to the slice shown as the gray rectangle. The current of course drops to 0 at the wire surface since we assume the wire is surrounded by an insulating medium. Conversely, in any planar slice of the wire which is perpendicular to the gray plane (and still parallel to the z axis), there are no eddy currents because any math loop in such a plane sees no flux. Suppose now that the field Bext has a positive linear gradient in the x direction, meaning that Bext has a larger amplitude toward the top of Fig *** and a smaller amplitude at the bottom. We can write, jω Bext(x,ω) = jω[ Bext0 + αx] => ext(x,t) = ext0 + x Write ext(x,t) = ext(t) + (t) x = ω ejωt [ Bext(0) + α(0) x ] We assume α(0) > 0 so the B field magnitude is larger at the top of Fig ** than at the bottom at t = 0. Below we shall assume a time such that ejωt = -1 so then both ext(t) and (t) are negative. Then the first equation above becomes, ∂xJz = σ ext = σ [ext0 + x ] => Jz(x) = σ [ext0 x + (1/2) x2 - (1/6) ] or Jz(x) = - σ [ | ext0| x + (1/2) | | x2 - (1/6) | | ] // for our time of interest where we have added a constant such that !Syntax Error, Idx Jz(x) = 0 for a wire of radius a = 1. Jeddy = Jz is now larger in the upper half of the gray rectangle than in the lower half, and we would expect then a pattern having this general shape, (e) Fig **. Eddy Currents in a quiet wire (I = 0) induced by an external B field with a gradient. The eddy currents are now larger on the side of the wire where the external field Bext is larger. In addition, one sees the eddy current in general to be larger near the surface of the wire and small in the interior. P.6 Eddy Currents induced in an current-carrying wire by an external B field Continuing with our series of toy examples, we now superpose a small amount of Figure (d) onto Figure (b) in order to represent the eddy current for a current-carrying wire which is placed in a small uniform external magnetic field Bext. The current I in the wire and the Bext are assumed to vary at rate ω. The result is something like this: (f) On the top there is some cancellation between the two eddy current patterns, while on the bottom there is reinforcement. Thus we arrive at another mechanism for the eddy current to be larger on one side of a wire than on the other side. Notice that Bext has no gradient in this example. P.7 Summary of Round Wire Examples 1. The eddy currents which a current-carrying round wire induces into itself vary with radius inside the wire, but are azimuthally symmetric. These eddy currents are interpreted as causing the skin effect. This situation is depicted in Fig b. 2. A quiet round wire in the presence of a spatially-uniform external B field will have induced eddy currents which are oppositely directed on the two sides of the wire (top and bottom side), but the absolute value of the current is symmetric on the two sides, as in Fig d. 3. If this quiet wire is placed in an external Bext field which has a gradient as shown above, then the absolute value of the current density will be larger on the side of the wire where Bext is larger, as shown in Fig e. 4. When a current-carrying round wire is placed in a uniform external Bext field, even though that field is uniform, the absolute value of the current density is larger on one side of the wire compared to the other side, as shown in Fig f. 5 When a current-carrying round wire is placed in a external Bext field which has a gradient, we again expect to have a top/bottom eddy current asymmetry which is a combination of the effects of items 3 and 4 above. The asymmetry will depend on the direction of the current in the wire and on the size and polarity of Bext. P.8 Eddy currents in Transmission Lines: The Proximity Effect Consider a transmission line composed of two round wires and we focus our attention on wire #1 as our Device Under Test. If wire #2 is far away, as in a wide-spaced twin-line, then Bext created by wire #2 is roughly uniform at the location of wire #1, and we then have the current asymmetry of Case 4 above. If wire #2 is close to wire #1, then Bext created by wire #2 will have a gradient over wire #1, and then we have the asymmetry combination Case 5 above where both effects must be considered. In the following drawing, we show in cross section two round wires both of which carry current I in the same direction, out of the plane of paper. The field Bext created by wire #2 is slightly stronger on the right side of wire #1 than on the left side, which of itself would argue for more eddy current on the right side of wire #1. However this effect is swamped by the Case 4 effect where we have cancellation of B fields on the right side of wire #1 and addition on the left side, so the B field is stronger on the left side of wire #1 and thus the current is larger there, as indicated by the lighter coloration. The current asymmetry increases as the two conductors get closer together because both the B field cancellation and addition are enhanced as Bext becomes larger are more comparable to the internal B field. If we imagine positive charge carriers coming out of the plane of paper in wire #1, the ones on the left side of wire #1 feel a Lorentz force a q v x B pushing them to the right, while those on right side of wire #1 feel an oppositely directed force pushing them to the left. But the charge carriers on the right are in a smaller B field and travel at a smaller velocity v since Jz = nqv is smaller. The net effect is that the charge carriers in wire #1 feel an overall force to the right and this force is transferred to the conductor ion lattice to maintain ρ = 0 causing the entire wire to be pushed to the right. The opposite happens inside wire #2 and the result is that wires with currents in the same direction attract each other. The fact that the current distribution in each wire is skewed away from the other wire is sometimes called the proximity effect, or current crowding. One effect of having a non-uniform Jz distribution is that the wires have resistance larger than their DC values. If the above two wires were two strands of a power transmission line cable carrying current in the same direction, the Ohmic loss in the strands is enhanced by this effect. The coloration patterns in our figure don't really illustrate the skin effect which is of course always present at any ω > 0, more strongly of course when δ < a. Even in power lines at 60 Hz where δ ~ 1 cm the skin effect does cause a waste of the wire interior since less current flows their. If a large round conductor is replaced by a set off smaller insulated round wires, this waste is reduced since the smaller wires each have more uniform current (see Litz wire). In a normal transmission line the currents are of course oppositely directed and the picture is different: Now the B field is larger on the right side of wire #1 so the eddy current is larger there causing the total current Jz to be larger there compared to the left side. The currents are now crowded on the side of each wire facing the other wire. The Lorentz force now causes the two wires to repel each other, reversing the argument given above. There is still Ohmic loss compared to DC since Jz is non-uniform. Closer wire spacing again results in increased asymmetry. Companies like "Monster Cables" advocate using their low-ohm expensive cables for driving audio speakers in order to offset the resistance increases due to both proximity and skin effects. If one ignores the proximity effect and assumes a uniform Jz, it is not hard to compute the total magnetic field B of the two conductors at any point (x,y) in the cross-section plane. In the following graphs we show |By| (red) as a function of x in the y = 0 slice through the two wires. The wires have radius 1/2 unit and center separation 3 units. Currents in same direction: Currents in opposite direction These graphs support the claims made earlier about where the total B field is large and small, and thus where the eddy current is large and small. Of course the graphs are approximate since Jz is in fact not uniform on each wire cross section. P.9 Quantitative Evaluation of Eddy Currents and The Proximity Effect Our presentation above has all been qualitative, and no method was given for computing the actual size of the eddy currents and thus of the proximity effect. G. Smith presents the following intriguing chart showing the strength of the effect versus conductor separation under unspecified conditions, We present in Chapter 7 a quantitative treatment of the skin and proximity effects which gives results similar to the above graph. P.10 Proximity Effect and Wire Resistance Consider a small differential volume rdθdrdz in one of the round wires (relative to a cylindrical coordinate system for that wire). Its cross sectional area is dA = rdθdr . This volume has resistance dR = "ρL/A" = ρ dz / dA and the current through this resistor will be dI = Jz(r,θ) dA . The Ohmic power generated in this tiny resistor is, from P = I2R, dP = (dI)2(dR) = [Jz(r,θ)dA]2 ρ dz / dA = ρ Jz(r,θ)2 dA dz . For the larger resistor consisting of length dz of the entire cross section we find then that P = ∫dP = dz ∫dA ρ Jz(r,θ)2 . The total current in the wire is I = ∫dA Jz(r,θ) and then from P = I2R the effective wire resistance of a cross sectional slice of wire of length dz is, R = = ρdz ∫dA [...] = !Syntax Error, Idr r !Syntax Error, Idθ [...] Adding some cancelling factors of A = πa2 we get, R = R/dz = [ ρ/A ] = Rdc = Rdc where Rdc is the DC resistance per unit length of the wire. Our notations <> and E() mean "expected value". In elementary probability theory one writes μx = E(X) // mean σx2 ≡ varx = E(X2) - E(X)2 = E(X2) - μx2 // variance; σx = standard deviation so that = = 1 + Thus, taking X = Jz we find this result for AC resistance per unit length, R = Rdc (1 + ) = 1 + = [1 + ] // loss At DC, Jz is constant across the cross section so its variance is 0 and the above says R = Rdc. For any other function Jz(r,θ) ≠ constant, one will have some variance σJz2 > 0 and then R > Rdc. Thus, the proximity effect increases the effective resistance of the wires in Fig XX, causing an increase in the Ohmic loss. Notice that the percentage proximity loss is independent of the current I. ******************* This reference should be added: G. Smith, "The Proximity Effect in Systems of Parallel Conductors and Electrically Small Multiturn Loop Antennas" (Harvard University Division of Engineering and Applied Physics, Technical Report No. 642, Dec 1971) see www.dtic.mil/dtic/tr/fulltext/u2/736984.pdf‎ . This document contains a long list references, see also the search engine at www.dtic.mil/dtic .