Confinement of high temperature plasmas

The production of useful energy by nuclear fusion has been the subject of extensive research in a number of different countries for approximately twenty years. A recent world-wide survey of various fusion research facilities and the fusion devices at those facilities is reported in "World Survey of Major Facilities in Controlled Fusion" (Amasa S. Bishop ed. 1970). However, although many different approaches have been suggested and tried, no one has yet succeeded in producing controlled nonexplosive energy in useful quantities by nuclear fusion.
Research on the generation of energy by nuclear fusion has centered on the problem of confining within a vacuum relatively low density gases, such as deuterium, which have been raised to temperatures sufficiently high for fusion reactions to take place. At the temperatures required for fusion reactions to take place, such gases become fully ionized and are referred to as plasmas, and the inability to confine such high temperature plasmas for a sufficiently long time to obtain more energy from fusion reactions within the plasma than is supplied to the plasma, has been the principal obstacle to the useful production of nonexplosive energy by nuclear fusion.
A number of arrangements have been proposed for confining high temperature plasmas within a vacuum, most of which are based upon using one or more magnetic fields which exert forces on the electrically charged particles composing the plasmas. Several of these magnetic confinement arrangements, sometimes referred to by the term "magnetic bottles," are disclosed in "Progress Toward Fusion Power," Fowler and Post, 215 Scientific American 27 (1966); "The Prospects of Fusion Power", Gough and Eastlund, 224 Scientific American 50 (1971); "The Road to Controlled Nuclear Fusion", Artsimovich, 239 Nature 18 (1972); and U.S. Pat. Nos. 3,009,080 granted to H. G. Loos on Nov. 14, 1961, and 3,088,894 granted to H. R. Koenig on May 7, 1963. However, none of these magnetic confinement schemes has been capable of maintaining the confinement of a plasma for more than brief periods which are too short to be of practical usefulness in the production of energy. The inability of known magnetic confinement schemes to maintain the confinement of plasmas for more than very brief periods has been attributed to the unstable nature of plasmas, and has led to the investigation of the possible instabilities of various plasma configurations, and of techniques for overcoming these instabilities. See, for example, Elementary Plasma Physics, Arzimovich (1963, 1965), Chapter 8.
One of the widely used techniques attempted for the confinement of plasmas is to form the plasma in a toroidal vessel so that the plasma is shaped into a ring shaped loop resembling a "doughnut". Examples of variations of this technique are the "Stellarator", "Scyllac" and "Tokomak" machines described in the above-mentioned articles, and the apparatus disclosed in the above-mentioned Loos and Koenig U.S. Pat. Nos. 3,009,080 and 3,088,894. Although this technique and its variations have been known and used for a number of years, long term plasma stability sufficient for the production of useful energy by nuclear fusion has not yet been achieved.
BRIEF SUMMARY OF THE INVENTION
This invention relates to the generation of energy by nuclear fusion, and in particular to the confinement of high temperature plasmas sufficiently long for the energy generated by fusion reactions within such a plasma to exceed the energy supplied to the plasma both to compensate for plasma energy losses, and, if necessary, to maintain the plasma at a sufficiently high temperature.
The present invention provides a solution to the problem of confining a high temperature plasma sufficiently long to obtain useful energy by nuclear fusion, based upon the discovery that there exists a class of toroidal magnetic confinement arrangements characterized by "structural stability". "Structural stability" is defined as that property of plasmas confined by the class of magnetic confinement arrangements of this invention, which makes such plasmas relatively insensitive, for sufficiently long periods of time, to approximations, tolerances and variations in the parameters of design, construction and operation of such arrangements.
The class of magnetic confinement arrangements of this invention is characterized by the provision of magnetic field lines at the surface of a toroidal plasma which satisfy the requirements of static structural stability. These requirements include the conditions that only a finite number of magnetic field lines are closed. In addition, in order to prevent plasma leaks, it is further required that the magnetic field be everywhere non-zero. Finally, for reasons explained in detail in the following description, it is required that the limit cycles of the magnetic field lines be even and non-zero in number and that the rotational transform angle of the magnetic field lines be identically and persistently zero. Magnetic field lines that satisfy these requirements may be achieved by a set of current carrying conductors properly arranged with respect to the exterior of the plasma so that the resulting field lines generated at the surface of the plasma by the currents within the conductors include an even finite number of closed toroidal field lines which divide the plasma surface into a corresponding even number of toroidal regions, where the poloidal field lines in each region are opposite in direction to the poloidal field lines on the adjacent region or regions and the rotational transform angle of all the magnetic field lines in all of the regions is zero. Toroidal plasmas confined by such magnetic field lines constitute a class of plasmas in which each of the members of the class is characterized by a non-convex poloidal cross-section, for example, a cross-section resembling either a "kidney bean", that is, an approximation of the Limacon of Pascal, or a "banana" shaped curve, and in which the poloidal cross-section of such plasma has an even finite number of stagnation points which correspond in number and location to the closed magnetic field lines on the surface of the plasma.

BRIEF DESCRIPTION OF THE DRAWINGS
The present invention will be fully understood from the following detailed description of illustrative embodiments thereof taken in connection with the appended drawings, in which like reference characters designate like elements throughout the several drawings:
FIG. 1 is a schematic perspective view partially in section illustrating the various current conductors and other elements embodying the principles of this invention;
FIG. 2A is a cross-section of a portion of the elements shown in FIG. 1 in which the non-convex cross-section of the plasma is necessary for static structural stability for the reasons given in the following detailed description;
FIG. 2B is a cross-section of a portion of elements alternative to those shown in FIG. 2A;
FIG. 3A is a perspective view showing in partial section the families of magnetic surfaces generated by the current conductors shown in FIG. 2A;
FIG. 3B is a pair of side views of a plasma confined in accordance with this invention, the upper of said pair of side views illustrating an outside view of the magnetic field lines on the outer or larger radius surface of said plasma, and the lower of said side views illustrating an inside view of the magnetic field lines and limit-cycles or closed magnetic field line orbits on the inner or smaller radius surface of said plasma;
FIG. 4A is a perspective view partially in section illustrating a toroidal plasma configuration in accordance with the principles of this invention;
FIG. 4B is a cross-section of the toroidal plasma shown in FIG. 4A;
FIG. 5 is a cross-section of an alternative toroidal plasma configuration in accordance with the principles of this invention in which the number of poloidal stagnation points is equal to four;
FIG. 6 is a cross-section of another toroidal plasma configuration in accordance with the principles of this invention;
FIGS. 7, 8 and 9 are graphs of assistance in explaining the theoretical principles in Part II hereof underlying this invention;
FIGS. 10 through 33 are graphs of assistance in following the descriptions in Part II, sub-parts A-M of the specification.
FIGS. 34, 35 and 36 illustrate a specific embodiment of the invention in which locations of wires and current values are shown.
To avoid unnecessarily complicating the drawings, the insulators or other means well known to those skilled in this art for maintaining the conductors in their proper positions have not been shown. For the same reason, power supplies shown in the drawings should be understood to be representative of a multiplicity of power supplies necessary to drive those several currents disclosed herein as required by the principles of this invention.

DETAILED DESCRIPTION OF THE INVENTION
PART I
APPARATUS EMBODYING THEORETICAL PRINCIPLES
Referring to FIGS. 1 and 2, a hollow current carrying conductor 10, which may be either a single integral conductor or a plurality of segmented conductors, defines a toroidal chamber 11 within which there is mounted a suitably non-conductive toroidal envelope 12 formed of an electrically insulating, magnetically permeable material such as quartz or porcelain. Both conductor 10 and envelope 12 have a common axis Z which will be referred to in this specification and in the claims as "toroidal axis" or "axis of the toroid". The interior of envelope 12 is divided by internal wall 120 into two gas-tight chambers 121 and 122, and chamber 121 is filled with a suitable gas such as deuterium at a relatively low pressure, for example, on the order of one-thousandth of atmospheric pressure. Within chamber 122 there are mounted toroidal or azimuthal current carrying conductors 123 and 126 also having axis Z as their common toroidal axis. Conductor 123 may take the form of either a plurality of toroidal current carrying wires arranged about a longitudinal axis L1 and also having axis Z as a common toroidal axis, as shown by circles enclosing a dot in FIGS. 1 and 2, or a single toroidal current carrying conductor centered on longitudinal axis L1 as shown in FIG. 3A. Conductor 126 is embedded in an insulator which also serves as a support for 123, and conductor 126 may take the form of a plurality of current carrying wires as shown in FIGS. 1, 2A and3A.
Mounted exteriorly of envelope 12 is a toroidal current carrying conductor 13 which may take the form of a plurality of toroidal current carrying wires illustrated by small circles each enclosing an "x" in FIGS. 1 and 2, all of which have axis Z as a common toroidal axis.
The gas within chamber 121 is first ionized to become a toroidal plasma, having Z as its toroidal axis, by any one of a number of well-known arrangements (not shown) such as those disclosed in "Controlled Thermonuclear Reactions", Glasstone & Lovberg, pp. 116-163 (1960). Following ionization of the gas within chamber 121 to form a plasma, current carrying conductors 10, 123 and 13 are respectively energized in the following manner by suitable capacitor banks 130, 131, 132 such as those shown in H. G. Loos U.S. Pat. No. 3,009,080, to confine the toroidal plasma within chamber 121 to have the configuration 14 illustrated in FIGS. 2A, 2B, 3A, 3B and 4A. Current is directed poloidally through conductor 10, that is, transversely around the toroidal chamber 11 enclosed by conductor 10 as shown by the arrows in FIGS. 2A and 2B. At the same time current is directed through conductors 123, 126 and 13, with current being directed through conductors 123 and 126 in a direction opposite to the direction of current through conductor 13. For example, current may be directed through conductor 13 in a counterclockwise direction around toroidal axis Z, as indicated in FIG. 1 by an "x" within the circles denoting the wires of conductor 13, whereas current is directed through conductors 123 and 126 in a clockwise direction around toroidal axis Z, as indicated in FIG. 1 by the dots within the circles denoting the wires of conductors 123 and 126.
As will be explained in detail in the following portions of this specification, the currents flowing through conductors 10, 123, 126 and 13 generate within the interior of envelope 12 sets or families of toroidal magnetic fields or surfaces 124 and 125, as indicated by the lines with arrows in FIG. 2A and by the examplary surfaces in FIG. 3A. The existence of two stagnation points, 4 and 5, in the poloidal component of the magnetic field is evident in FIG. 3A from the directions of the arrows which indicate the direction of the magnetic field lines produced by the currents in the plasma surface and in the wires 123 and 126. These currents also induce magnetic field lines 19 on the surface of plasma 14, as shown in FIG. 3B, the magnetic surfaces 124, 125 and the magnetic field lines 19 having certain properties required by this invention as defined below.
The magnetic field contributed by the current flowing through conductor 10 tends to displace the magnetic field lines 19 on the surface of plasma 14 inwardly toward its longitudinal axis L2, thereby tending to compress the plasma about its longitudinal axis L2, as well as to curve the magnetic field lines 19 on the surface of plasma 14 about toroidal axis Z. Currents in the wires 10, 123, 126 and 13 are adjusted pursuant to the criteria of static structural stability explained herein, thereby to yield the final equilibrium configuration of the magnetic field.
The magnetic field contributed by the current flowing through conductors 123 and 126 adjusted in accordance with the principles of this invention serves to provide a region on the inner surface of plasma 14, that is, on the surface nearer toroidal axis Z and the smaller radius or internal wall 120 of envelope 12, on which magnetic field lines 19 are curved or pushed toward the longitudinal axis L2 of plasma 14 to form a "dimple" so that the poloidal cross section of plasma 14 is non-convex. This region is illustrated in the lower side view of FIG. 3B by the area between lines 15 and 16. The non-convex poloidal cross section may take any one of several forms, for example, a "kidney bean" shape 141 as shown in FIG. 2A, 2B, 3A and 4B, a modified "kidney bean" shape 142 as shown in FIG. 5, or a "banana" shape 143 as shown in FIG. 6.
The magnetic field generated by conductor 13 acts upon plasma 14 to prevent outward movement of plasma 14 away from torodial axis Z toward the outer wall or larger radius side of chamber 121.
The combined effect of the magnetic field provided by currents through conductors 10, 123, 126, and 13 serves to confine plasma 14 so as to satisfy the requirements of static "structural stability" explained in detail in Part II of this specification.
Referring now to FIG. 3A, the two families or sets of magnetic fields or surfaces 124, 125 which are generated within envelope 12 by current through conductors 10, 123, 126 and 13, are characterized by the following properties. The family of magnetic fields 124 is nested toroidally about conductor 123, and the other family of magnetic fields 125 is nested toroidally between envelope 12 and the outer boundaries of both the family of magnetic fields 124 and plasma 14. The magnetic field lines enter and leave the vacuum volume surrounding plasma 14 from the conductor 126. There is discontinuity in the tangential component of the magnetic field at the point where the magnetic field lines enter and leave from conductor 126, and the magnitude of this discontinuity is determined by the specific values of the currents in the conductors 10, 123, 126 and 13. It will be observed in FIGS. 3A and 3B that the poloidal component of the magnetic field at the surface of plasma 14 is not uniform in direction, but is characterized by two oppositely directed magnetic fields. Thus on the inner segment 465 which comprises arc 4 to 6 plus arc 6 to 5 of the boundary .GAMMA. of the poloidal cross section of plasma 14, that is, on the boundary of the poloidal cross section of plasma 14 which is nearer conductor 123 and the toroidal axis Z of plasma 14, the poloidal mangetic field line has a direction which is opposite to the direction of the poloidal field line on the outer segment 475 which comprises arc 4 to 7 plus arc 7 to 5 of .GAMMA. which is adjacent the outer wall of envelope 12 and further away from the toroidal axis Z of plasma 14. As a result, the poloidal components of the mangetic field lines 19 on the surface of plasma 14 are divided into two groups, each following the direction of the next adjacent one of the two sets of magnetic fields 124, 125 generated within chamber 121. As shown in FIGS. 3A and 3B, the poloidal components of the magnetic field lines 19 on the inner segment 465 of boundary .GAMMA. have the same direction as adjacent magnetic surface family 124, while the poloidal components of the magnetic field lines 19 on the outer segment 475 of boundary .GAMMA. follow the direction of the adjacent magnetic surface family 125.
As explained in detail in Part II of this specification the rotational transform angle of the magnetic field lines 19 on the surface of plasma 14 is required by the principles of this invention to be identically and persistently zero for all sufficiently small deviations in the parameters of the overall system. This requirement is illustrated in FIGS. 2A and 3B, in which the magnetic field lines 19 on the plasma surface are divided into two regions by magnetic field lines 15 and 16 on the surface of plasma 14. Field lines 15 and 16 are "closed", that is, they describe a closed circular, toroidal orbit on the surface of plasma 14. In accordance with principles of this invention as set forth in Part II of this specification, the number, N, of closed magnetic field lines such as 15, 16 must be even, that is, N=2, 4, 6, . . . , and these closed field lines divide the surface of plasma 14 into a corresponding number of N toroidal regions. Thus as illustrated in FIG. 3B (inside view), the two closed magnetic field lines 15 and 16 divide the surface of plasma 14 into a corresponding pair of toroidal regions. It is explained in Part II how to compute the current strengths in conductors 10, 123, 126 and 13 in order to achieve a magnetic field on the surface of plasma 14 with these required properties.
Within the plasma surface regions between the closed magnetic field lines, all other magnetic field lines describe non-closed toroidal paths which asymptotically approach or depart from the closed field lines. This is illustrated, for example, in FIG. 3B (inside view), in which the non-closed magnetic field lines 19 asymptotically depart from closed field line 15 and asymptotically approach closed field line 16 but do not cross either of closed lines 15 and 16. FIGS. 2A and 3B also illustrate that the path of a non-closed field line 19 is displaced by a finite angle, .delta., during each orbit around the circumference of plasma 14. This angle cannot exceed 2.pi. or 360.degree. during any given orbit around the circumference of plasma 14 because no non-closed magnetic field line 19 can cross closed lines 15 and 16. The total displacement angle .delta..sub.n after n circuits of a magnetic field line 19 around the circumference of plasma 14 cannot exceed 2.pi. for the same reason.
By definition, the rotational transform angle of magnetic field lines 19 is the limit of the ratio of total displacement angle .delta..sub.n to 2.pi.n, for an arbitrarily large value of n, that is ##EQU1## Because .delta..sub.n cannot exceed 2.pi. in order to comply with the structural stability principles of this invention, the limit ratio .delta..sub.n /2.pi.n, for n arbitrarily large, is identically zero and remains persistently zero for non-closed magnetic field lines 19. In addition, because all of the closed field lines have zero rotational transform angles by virtue of their closed nature, the entire class of both closed and non-closed magnetic field lines on the surface of plasma 14 is characterized by a zero rotational transform angle. The requirement of a finite even number of closed field lines to which all other field lines on the surface are asymptotic is an important feature of the present invention. Hyperbolic stagnation points of the poloidal field must be produced at a critical radius disclosed herein by a specific configuration of external currents computed in accordance with the principles disclosed herein.
It is to be understood that a plasma characterized by structural stability, as defined in this specification, may be obtained by means for generating magnetic fields other than those specifically disclosed herein, but which follow the principles of this invention. Thus there may be positioned within the vacuum volume surrounding the plasma, conductors which are spatially arranged in locations other than or in place of all or a portion of conductors 13, 123, and 126, for example, as shown in FIG. 2B, a suitably arranged plurality of conductors 21, 22 and conductor pairs 23a, 24b which provide magnetic dipoles in the poloidal field. Either well-known analytical or numerical techniques may be employed to select and locate appropriate conductors; for example, see Part II of this specification and sub-parts C and D thereof.
It is also to be understood that although the arrangements shown in FIGS. 1 and 2A disclose conductors at or coincident with flux surfaces or flux singularities, for example, the outer member of the family of surfaces 125 or the inner member of 124, surfaces other than flux surfaces may be selected, in which case appropriately chosen conventional magnetic dipole devices may be employed as shown in FIG. 2B to provide the magnetic field which would otherwise be generated by conductors at the flux surfaces.
Example
The poloidal cross-section of a toroidal plasma provided with topological stability in accordance with the principles of this invention is illustrated in FIGS. 34, 35, and 36. The plasma has an outside diameter of 1.02 meters and an overall height of 0.045 meters. The boundary, .GAMMA., of the plasma cross-section is described by the set of position values given in Table I below. These position values were obtained in accordance with the principles of this invention by locating the conductors 13 to produce a tangent magnetic field constant in magnitude everywhere on the surface of the plasma.
The plasma equilibrium configuration for the embodiment shown in FIG. 35 may be obtained by magnetic induction, in which case the magnetic pressure is composed of two contributions, one toroidal and the other poloidal. ##EQU2## The toroidal contribution is: ##EQU3## where B.sub.o and r.sub.c are initial values of the total equilibrium magnetic field magnitude on the plasma surface and the minimum radius of the plasma equilibrium boundary cross-section, respectively. The poloidal contribution is: ##EQU4## where the surface current density K is that of the boundary condition
n.times.(B.sub.+ -B.sub.-)=.mu.K
(see Stratton, Electromagnetic Theory, McGraw-Hill, 1941.) The surface current density K.sub.n corresponding to the position of the n.sup.th conductor on the plasma surface is ##EQU5## where J.sub.n is the current induced in that conductor and .DELTA.l.sub.n the length of boundary corresponding to that conductor allotted to it.
The locations and currents of the external conductors for the embodiment shown in FIG. 35 are given in Table II. FIG. 36 illustrates the position of the equilibrium plasma boundary r relative to the locations of the conductors 13a, 13b+, 13b-, 13c+, 13c-, 13d+, 13d-, 13e specified in Table II. (To avoid unnecessarily complicating the drawings of FIGS. 34 and 36, the insulators or other well known means for maintaining these conductors to their proper positions have not been shown). These currents correspond to a total magnetic field at the plasma of 200,000 Gauss. Lower magnetic fields would require proportionally lower currents.
TABLE I______________________________________Point PositionsOf Plasma Boundary .GAMMA.Shown in FIG. 35r (meter) z (meter)______________________________________.5251 .+-..0007.5250 .+-..0023.5248 .+-..0039.5245 .+-..0054.5240 .+-..0070.5235 .+-..0084.5229 .+-..0099.5222 .+-..0113.5214 .+-..0127.5205 .+-..0140.5195 .+-..0152.5184 .+-..0164.5173 .+-..0175.5161 .+-..0185.5149 .+-..0195.5136 .+-..0204.5122 .+-..0213.5108 .+-..0220.5094 .+-..0227.5080 .+-..0233.5065 .+-..0238.5049 .+-..0243.5034 .+-..0246.5018 .+-..0249.5003 .+-..0251.4987 .+-..0252.4971 .+-..0252.4955 .+-..0251.4940 .+-..0249.4924 .+-..0247.4908 .+-..0244.4893 .+-..0240.4878 .+-..0236.4863 .+-..0232.4847 .+-..0227.4832 .+-..0225.4832 .+-..0209.4840 .+-..0195.4847 .+-..0181.4852 .+-..0166.4857 .+-..0151.4862 .+-..0136.4866 .+-..0120.4869 .+-..0105.4872 .+-..0089.4874 .+-..0074.4876 .+-..0058.4878 .+-..0042.4879 .+-..0027.4878 .+-..0011______________________________________
TABLE II______________________________________Locations And Currents OfExternal Conductors Shown in FIG. 36Wire r (meters) z (meters) current (amps)______________________________________13a .45 0 +600,00013b+ .45 +.04 -130,00013b- .45 -.04 -130,00013c+ .47 +.07 -130,00013c- .47 -.07 -130,00013d+ .54 +.06 -240,00013d- .54 -.06 -240,00013e .58 0 -240,000______________________________________
PART II
THEORETICAL PRINCIPLES
The following nomenclature will be used in this portion of the specification and sub-parts A-M thereof:
______________________________________NOMENCLATURE______________________________________.fwdarw.B magnetic field vector.fwdarw.r position radius vector.fwdarw.j electric current vectorB magnetic field strength (length of .fwdarw.B or .fwdarw.B.multidot. .fwdarw.B).mu. magnetic permeabilityB.sup.2 /2.mu. magnetic pressurep plasma fluid pressure.beta. ratio of plasma fluid pressure to external magnetic pressure.SIGMA. surface (compact orientable 2-manifold) forming interface between plasma and vacuum, i.e. plasma boundary.fwdarw.n normal vector to .SIGMA.s arc length along field line.alpha. pertaining to positive limit set of a field line.omega. pertaining to negative limit set of a field line(x.sub.1, x.sub.2) local coordinates on .SIGMA.(X.sub.1, X.sub.2) components of .fwdarw.B on .SIGMA. local coordinate system.gamma. closed field line on .SIGMA. (topological circle)(r, .theta., z) cylindrical coordinates.psi. magnetic flux function.GAMMA. intersection of .SIGMA. with plane .theta. = O.differential. partial differential operatord differential operator.gradient. gradient operator.DELTA. Laplace's operator (.DELTA. = .gradient. .multidot. .gradient.).DELTA.* Stokes' operatorV volumeS surface arear.sub.c critical radius.epsilon. indicator______________________________________
In the following description of the theoretical principles underlying this invention, it will be understood that although there are two distinct aspects of structural stability of toroidal plasmas, static and dynamic, only static structural stability will be considered in detail. For a toroidal plasma characterized by static structural stability in accordance with the principles of this invention, dynamic structural stability can be achieved by well-known techniques such as systematic theoretical multi-channel design procedure of the type used in the design and development of aerospace-vehicle stabilization systems. Accordingly, unless expressly stated otherwise, the term "structural stability" in the following description will refer only to static structural stability.
Consider a plasma confined in a certain finite bounded region of space by a magnetic filed B, where the plasma contains no internal conductors of current and where there is litle or no magnetic field inside the plasma. At the plasma boundary or surface .SIGMA. forming the interface of the plasma and the magnetic field B, it is required that the total pressure p.sub.t, which is the sum of the plasma fluid pressure, p, and the magnetic pressure, B.sup.2 /2.mu., where .mu. is the permeability, that is,
p.sub.t =p+B.sup.2 /2.mu.,
must be continuous across the plasma-vacuum boundary. That is, the total pressure must be the same on each side of surface .SIGMA. (see Part II, sub-part A thereof). Using the subscript "i" to denote the pressure inside the plasma and the subscript "o" to denote pressure outside the plasma, one has ##EQU6## Because the plasma fluid pressure p.sub.o is zero outside of the plasma, equation (2) can be written The quantity .beta. is defined as the ratio of plasma fluid pressure inside the plasma, p.sub.i, to the magnetic pressure outside the plasma, B.sub.o.sup.2/ 2.mu., that is,
.beta.=p.sub.i /(B.sub.o.sup.2 /2.mu. ). (4)
It is evident from its definition that the universal parameter .beta. lies between zero and one. By requiring that there be little or no magnetic field inside the plasma (which is achieved for a high conductivity plasma, by induced currents in the plasma surface which prevent the penetration of magnetic field to the interior; this may be achieved by rapidly turning on the currents in the external conductors), equation (3) can be rewritten ##EQU7## and by substitution of the result of equation (5) in equation (4) it is evident that .beta. will approach unity under these conditions. Consideration of transport processes, particularly finite resistivity, shows that the requirement of dynamic stability of the plasma makes it preferable, but not necessary, that the value of .beta. approach unity (high .beta. values) as closely as possible. It is to be understood that the principles of this invention are equally applicable to low and medium values of .beta.. With unity .beta. under equilibrium conditions the fluid pressure on the interior of the plasma will have the same non-zero value everywhere, so that there will be a general steady fluid flow within the plasma (see Part II, sub-part H). Then according to equation (5), the magnetic pressure outside the plasma is non-zero everywhere on the boundary or surface of the plasma.
It is also required that the plasma surface .SIGMA. have no cusps or edges, and that the surface must be closed, that is, .SIGMA. cannot have any topological boundaries. It is further required that .SIGMA. be orientable in the sense that it has a unique inside and outside, that is, no part of .SIGMA. can interpenetrate another part of .SIGMA., or in mathematical terminology, that .SIGMA. be a compact orientable 2-manifold.
Another boundary condition for .SIGMA. (see Part II, sub-part A) is that the normal component of the magnetic field outside the plasma be continuous across .SIGMA., which implies that the magnetic field at the surface of the plasma is tangent to .SIGMA. and therefore forms a non-vanishing tangent magnetic vector field everywhere on .SIGMA..
It can be shown by the Lefschetz-Hopf fixed point theorem that by requiring .SIGMA. to be a compact orientable 2-manifold without any fixed point and the normal component of the magnetic field outside the plasma to be continuous across .SIGMA., then .SIGMA. must be topologically a torus, that is, the topological product of two topological circles.
The next requirement is that the confinement of a toroidal plasma be statically "structurally stable", that is, such a plasma must have equilibrium characteristics that are insensitive to small variations in those system parameters which determine the qualitative or topological properties of the magnetic field line configuration or portrait on .SIGMA. so that this portrait is unchanged by such variations. The necessary and sufficient criteria for structural stability of an arbitrary vector field on a topological torus have been set forth by M. Peixoto (See references in Foundations of Mechanics, Abraham (1967) and by V. A. Pliss (as expressed by the terms "coarseness" in Nonlocal Problems in the Theory of Oscillations (1964, 1966); see also, Systemes Grossiers", Andronov & Pontriagin, 14 Dokl. Akad. Nauk. SSSR (1937):
(a) there is only a finite number of singularities, all generic;
(b) The .alpha.- and .omega.-limit sets of every trajectory can only be singularities or closed orbits;
(c) no trajectory connects saddle points;
(d) there is only a finite number of closed orbits, all simple.
"Singularities" in the above criteria (a) and (b) refer to points where the vector field vanishes, and such points are usually called "stagnation points" in the art to which the present invention relates. A "saddle point" is a type of singularity also referred to as a "hyperbolic stagnation point".
Applying these criteria of structural stability to a plasma in a magnetic field under the conditions set forth above, in which the plasma boundary corresponds to the topological torus and the magnetic field corresponds to the vector field, criteria (a) and (c) are satisfied because the magnetic field is non-zero everywhere on the plasma surface .SIGMA. which means that there are no singularities.
"Trajectories" and "orbits" refer to the arc-length parameterization of the magnetic field lines obtained by integrating on .SIGMA. the ordinary nonlinear vector differential equation ##EQU8## The ".alpha." or ".omega." limit set of these magnetic field lines is the totality of limit points of the lines in the limit of the arc-length parameter s respectively approaching -.infin. or +.infin.. Because there are no singularities, by criterion (b) the limit sets must be closed or endless orbits, that is, the limit sets must be limit-cycles or closed orbits which are in the .alpha.- or .omega.- limit sets of other distinct magnetic field lines.
Criterion (d) requires that the number of closed orbits be both finite and "simple". In order to be "simple", ##EQU9## where the integrand of equation (6) is the two-dimensional divergence of the magnetic field on the plasma surface .SIGMA., in which x.sub.1 and x.sub.2 are local coordinates on the plasma surface, X.sub.1 and X.sub.2 are the corresponding components of the magnetic field line at those coordinates, and .gamma. is the closed orbit under consideration.
It is shown in Part II, sub-parts E, F, and G below that by combining the criteria of structural stability with Maxwell's Equations and Bernoulli's Theorem for steady fluid flow, as exemplified by the considerations underlying equation (3) above, the magnetic field line configuration on the surface of the plasma must satisfy the following criteria in order to achieve structural stability:
(a) the rotational transform angle of the magnetic field lines on the surface of the plasma is identically zero;
(b) the number N of closed orbits is finite and even (N=2, 4, 6, . . . ); and
(c) the closed orbits are not only simple but also must occur in pairs in which the members of each pair are asymptotically stable for increasing or decreasing arc-lengths (once a particular assignment of the sence of increasing arc length has been made).
Application of these criteria to the simplest case of two closed magnetic field lines (which is the non-redundant case because all of the other cases for N even require the duplication of certain elements), makes it apparent that one of the two closed field lines is asymptotically stable, that is, neighboring field lines approach it asymptotically, and the other closed field line is asymptotically unstable in the sense that neighboring field lines move away from it asymptotically. ("Asymptotically" as used here refers to parameterization by arc-length, not by time.) See, for example, closed lines 15 and 16 and the neighboring field lines 19 in FIG. 3B. As shown in Part II, sub-part G below, the requirement that the two field lines with closed orbits be simple, is equivalent to the requirement that the "flow" induced on the plasma surface .SIGMA. by arc-length parameterization of the field lines not be, on the average, area-preserving in the immediate neighborhoods of the two closed field lines.
It will now be shown that under the requirements set forth above, the equations for equilibrium is magneto-fluid dynamics have a unique type of solution in the form of a class of physically realizable toroidal plasma configurations. This class of physically realizable toroidal plasma configurations is characterized by the following properties:
(1) The two closed orbit magnetic field lines arc two limit cycles on the plasma surface, one asymptotically stable and the other asymptotically unstable.
(2) The poloidal cross-section of the plasma is characterized by two stagnation points in the poloidal magnetic field on the plasma surface, each stagnation point corresponding to the poloidal cross-section of one of the two limit cycles on the plasma surface.
(3) The poloidal cross-section is non-convex.
No plasma confinement process hitherto advocated has this characteristic of static structural stability including non-area-preserving flow in the neighborhood of each (isolated) closed field line, since none uses the possiblity of asymptotic field lines on the plasma surface, and in particular, in prior plasma confinement devices no closed field line is isolated. This does not contradict the volume-preserving nature of the corresponding vacuum field flow, since the nearest vacuum magnetic shells are neither nested with, nor, in the vacuum, homologous to .SIGMA.. As explained previously in connection with FIG. 3A, it has been discovered that there is a class of physically realizeable arrangements for magnetic confinement of plasmas within a vacuum, consisting of at least two one-parameter families of nested toroidal magnetic surfaces each immediately adjacent to at least a portion of the plasma and a third single member magnetic surface family constituting the plasma-vacuum interface, where none of the first two families of surfaces is homologous to the third single member magnetic surface family and where the third single member magnetic surface family is simultaneously characterized by an identically zero rotational transform angle. A specific set of families of nested magnetic surfaces satisfying these criteria is illustrated in FIG. 3A, in which family 124 is homologous to conductor 123 (disregarding conductor 126), and family 125 is homologous both to envelope 12 and the geometrical union of family 124 with the surface of plasma 14.
The terms "homologous", "union" and "homotopic" used in this specification have the respective meanings well-known in the discipline of mathematical topology. In addition, it will be understood that as used throughout this specificaton and in the claims, the terms "vacuum" or "vacuum volume" mean the geometrical volume surrounding a plasma, regardless of the presence or absence of magnetically permeable matter and exclusive of electrical conductors, within which there is a magnetic field which confines the plasma to be structurally stable in accordance with the principles of this invention.
The following properties of the practical realization of this invention in accordance with the theoretical principles of this inventin are all natural and reasonable from the viewpoint of practical engineering:
(a) finite plasma volume;
(b) connectedness of the plasma volume;
(c) local homotopic triviality of the plasma volume (in particular, no internal conductors);
(d) closure of the device;
(e) plasma surface is a compact orientable manifold;
(f) non-redundancy of limit cycles;
(g) axial symmetry;
(h) static structural stability.
Under these requirements the standard equations of equilibrium have a unique type of solution, radically different from those hitherto given serious consideration in this art. A mathematicl framework suitable for specification of this solution is presented in detail in Part II, sub-part B. The following are the essential elements of the mathematical framework for this type of solution.
Let z be a coordinate along the toroidal axis of symmetry, and let (r, .theta.) be polar coordinates in any plane through the z-axis (that is, any poloidal plane). Then (r, .theta., z) constitutes a standard system of cylindrical coordinates, and the vector B in this coordinate system has components (B.sub.r, B.sub.o, B.sub.z). Without loss of generality, the vacuum field may be characterized by arbitrary constants c and p.sub.o and a scalar magnetic flux function .psi.=.psi.(r, x) such that ##EQU10## The function .psi. is well-known to those versed in the present art, as is the fact that its gradient, .gradient..psi., is necessarily normal to the surface .SIGMA.. The constant c is determined by the current in the external poloidal conductors (or, equivalently, along the z-axis). The constant p.sub.o denotes the equilibrium plasma fluid pressure.
The constant c completely determines the toroidal magnetic field. The poloidal field is uniquely determined by the constant p.sub.o and the function .psi.; but, as will now be shown, .psi. is uniquely determined by choice of c, p.sub.o, and the geometry of .SIGMA..Hence, finally, once c and p.sub.o have been specified, the vacuum field is uniquely determined solely by the choice of the geometry of .SIGMA..
This situation is illustrated in FIGS. 4A and 4B where 14 denotes the topological torus representing the plasma volume, and .SIGMA. denotes the surface of the plasma volume. The poloidal cross-section of this torus is 141 and .GAMMA. denotes the poloidal cross-section of just the surface .SIGMA. and is the curve
.GAMMA.: .psi.(r, z)=0. (8)
On .GAMMA. the pressure balance condition (see Part II, sub-part B) is:
.vertline..gradient..psi..vertline..sup.2= 2.mu.p.sub.o r.sup.2- c.sup.2. (9)
Clearly, for small enough r, the right-hand side of equation (9) is <0, whereas the left-hand side is always .gtoreq.0. Therefore, there is a minimum r, r.sub.c, at which equation (9) can be satisfied, that is, there is a minimum r, r.sub.c, at which a physically real solution exists. Thus, as shown in FIGS. 4A, 4B and 5, the plasma volume must always lie to the left of the vertical line r=r.sub.c. The critical radius r.sub.c is given by ##EQU11## and the boundary condition which, together with equation (8), determines .psi. on the exterior of .GAMMA. as a solution of Stokes' Equation, is ##EQU12##
If .GAMMA. lies wholly to the left of r.sub.c, then a particular value of the .+-. sign choice in equation (11) must be made, for .differential..psi./.differential.n cannot then vanish anywhere on .GAMMA.. But this implies that the vector field on .SIGMA. necessarily has a non-zero rotational transform angle, which violates requirement (a) of the structural stability criteria of the present invention.
By way of example, if, as in FIG. 9, one were to choose as .GAMMA. the boundary C of an arbitrary convex body which touches r.sub.c at just one point, 18, then by equation (11), the poloidal field must have a stagnation point at that point because there is a stagnation point at any point where .GAMMA. touches r.sub.c. The point 18 must be an elliptic singularity, implying the presence of a solid conductor at this point, and if the plasma should touch an external conductor at point 18, this would ruin the thermal insulation required for controlled fusion reactions. However, by producing oppositely directed sheet-currents on the left and right hand sides of .GAMMA. and by allowing the square-root in equation (11) above to change sign at a point where .GAMMA. touches r.sub.c, the requirement of zero rotational transform is satisfied by converting the singularity from an elliptic one to a hyperbolic one; it then becomes evident that a second point of contact with r.sub.c, which implies a non-convex poloidal cross-section, will permit the sign to be changed back again to its original value, which is necessary, by equation (11), in order that .differential..psi./.differential.n shall have no discontinuities on .GAMMA.. (Otherwise Stokes' Equation would have not have a physically realizable solution). The resulting "kidney-bean" or non-convex poloidal cross section 141 shown in FIGS. 4A, 4B, together with its associated magnetic field, is the consequence of not only the choice of a non-convex poloidal plasma cross-section 141 but also the choice of oppositely directed toroidal sheet-currents which permit .differential.104 /.differential.n to change sign as indicated at points 4 and 5 in FIG. 4B.
Thus by choosing .GAMMA. to have a non-convex shapethat touches the line r=r.sub.c at two points, 4 and 5, the requirement of two limit-cycles on .SIGMA. is achieved. Since, at 4 and 5, r=r.sub.c, then by equation (11), .vertline..gradient..psi..vertline.=0, and this gives B.sub.r= B.sub.z= 0. Thus the only magnetic field component present at 4 and 5 is toroidal (azimuthal), and so, as shown in FIG. 3B, the two field lines 15 and 16 respectively passing through points 4 and 5 are circles on .SIGMA., parallel to the r.theta. -plane, with centers on the z-axis. On the segment 465 of .GAMMA. we can take .differential..psi./.differential.n>0 and on segment 475, .differential..psi./.differential.n>0 (or vice versa).
The postulated normal discontinuity in B across .SIGMA. is achieved by surface current on .SIGMA.. These currents are oppositely directed in the two regions on .SIGMA. defined by the two limit cycles indicated in FIG. 3B by lines 15 and 16. In the redundant case of N>2 limit cycles, there are N/2 pairs of closed orbit magnetic field lines which divide the plasma surface into N adjacent toroidal regions or N/2 pairs of adjacent toroidal regions. Accordingly, the non-convex poloidal cross-section of such a surface has N poloidal stagnation points corresponding in number and geometrical location to the N limit cycles; for example, FIG. 5 illustrates N=4 stagnation points. Similarly, there are N toroidal sheet currents, one on each of the N adjacent toroidal regions, where each toroidal sheet current is opposite in direction to the toroidal sheet current on the two immediately adjacent toroidal regions.
The curve 20 in FIG. 4B is a "flux surface", i.e. a surface .psi.=constant. Toroidal currents appropriately placed on such a surface, together with induced toroidal plasma surface currents, generate the vacuum poloidal field; additional, oppositely directed toroidal currents are needed at each of points 8 (It will be shown below that these points lie along a branch cut.). However, one could, if desired, replace the current at points 8 by a family of toroidal currents appropriately placed on one of the nearby flux surfaces surrounding the outermost one of the points 8, such as 9. Outside the surface generated by axial rotation of curve 20, strong poloidal currents (not shown), together with induced poloidal surface currents, produce the toroidal magnetic field in the vacuum.
Existence of an analytic solution to Stokes' Equation in the neighborhood of .GAMMA. is proved by the Cauchy-Kovalevski Theorem, provided that .GAMMA. is an analytic curve, and provided that .differential..psi./.differential.n is analytic on .GAMMA.; this last fact, underlies the physical realizability of the principles of this invention, which is proved in Part II, sub-part B.
For simplicity of calculation one may take .GAMMA. to be a Limacon of Pascal, although this may not lead to an optimal configuration because the resulting branch cut 8 turns out to be uncomfortably near to the plasma. Another possibility is to choose .GAMMA. to be a "banana" shaped curve as shown in FIG. 6, which can be done utilizing cap-cyclide coordinates as shown in Part II, sub-part I. In this case one finds two possible singularities, as at 61 and 62. However, from the point of view of the present objective, namely plasma confinement, there is no essential difference between the choice of .GAMMA. disclosed in FIGS. 4B and 6.
In order to verify the "simplicity" of the limit-cycles through 4 and 5 as shown in FIG. 3B (as required for structural stability), one can represent the solution for the magnetic field in a neighborhood of either limit-cycle by a power series expansion of B and .GAMMA. convergent near the limit-cycle as proved in Part II, sub-part G.
The above are the principal considerations in showing the existence and uniqueness of a solution satisfying the last-named requirements (a, b, c, d, e, f, g, h). If the condition (f) of non-redundancy is relaxed, then solutions of the type shown in FIG. 5 are possible where the number of limit cycles and therefore stagnation points is increased by a multiple of two.
Note that the vacuum volume is at least triply-connected, by which we mean that there exists in said volume at least three independent non-homologous, geometrical cycles, including one which is homotopic to a point. These cycles are, for example, 41 in FIG. 4A, .GAMMA. and 42 in FIG. 4B. Thus conditions (a, b, c, d, e) and (g, h) imply a plasma configuration which is a topological torus, axially symmetric with a non-convex poloidal cross-section and with the associated vacuum volume at least triply-connected (unless dipoles are used shown in FIG. 2B).
Using established terminology of this art, in which a "Tokomak" is referred to as low-beta toroidal plasma monopole, and a "Gulf Energy and Environmental Systems Doublet" is referred to as a low-beat toroidal plasma monopole doublet (see "The Tokomak Approach in Fusion Research", Coppi and Rem, 227 Scientific American 65 (1972)), the preferred form of the present invention may be characterized as a high-beta toroidal plasma dipole.
Turning now to the subject of dynamic structural stability, or equivalently, of (dynamic) asymptotic stability of the plasma pseudo equilibrium (cf. Part II, sub-part J), all five subjects of electrical resistivity, heat conduction, compressibility, viscosity, and the Maxwell "displacement current" (effect of finite rather than infinite speed of electromagnetic waves) are to be considered simultaneously, rather than one or a few at a time, which is all that has been hitherto accomplished in the art. The reduction of the eleven simultaneous partial differential equations of EMTFD (electromagnetothermofluiddynamical) wave motion to a simple scalar static equation, namely Helmholtz's equation for the geometrical shapes of the standing modes of sound waves, plus a seventh-degree polynomial dispersion relation for the complex frequencies of the plasma's normal modes of vibration, renders the physical nature of the plasma's macroscopic oscillations quite understandable. (A special case of the above analysis may be found by taking the determinant of simultaneous homogeneous linear equation (2.60) and (2.61) of R. J. Tayler's "Plasma Waves And Oscillations", p. 113 of Physics of Hot Plasmas, edited by Rye and Taylor, Oliver and Boyd, 1970.). There will be inside the plasma principally four kinds of three-dimensional waves, namely
(i) entropy waves;
(ii) "slow" magnetosonic waves;
(iii) "fast" magnetosonic waves;
(iv) electromagnetosonic waves.
Waves (ii) and (iii) are modified sound waves and modified Alfven waves, or vice versa (depending on which is slower), while waves (iv) are modified electromagnetic waves. One novel conclusion of this analysis is that, among hitherto considered tubular plasma column-confinement methods, only those in the category of "infinitely long" (that is, open-ended) theta pinches (that is, of zero rotational transform) possess dynamic structural stability. Another result (see FIGS. 7 and 8 and Part II, sub-part K) is that it is easier to achieve dynamic structural stability for high beta (near .beta.=1) than for low beta (near .beta.=0). In particular, the analysis below and the graphs in FIGS. 7 and 8 indicate that over-idealized "Energy Principle" calculations, in which four of the above-mentioned physical effects are ignored, yield physically reliable results only regarding dynamic instability at high beta. Thus Teller's celebrated geometrical criterion for the dynamic instability of diamagnetic plasmas may be employed. Use of this criterion in the manner of Pfirsch and Wobig (see Pfirsch and Wobig, Proceedings, Culham Conference on Nuclear Fusion, Vol. I, 1966, p. 757) shows that interchange or global flute instabilities cannot occur because in the appropriate sense of "minimum average B" the present configuration is an optimal toroidal "magnetic well". Obviously the outer region of the plasma surface (between points 43 and 44 in FIG. 4B) has unfavorable curvature and may, for example, suffer local ripple instabilities of the type known in this art as "local flute instabilities" (see Pfirsch and Wobig, Proceedings, Culham Conference on Nuclear Fusion, Vol. I, 1966, p. 757). Such instabilities grow until the assumptions of linearizeability made in their analysis are no longer valid, but need not grow indefinitely. The unlimited growth of such instabilities can be suppressed effectively in at least two independent ways:
(i) by choosing a large enough toroidal aspect-ratio, (i.e. ratio of major to minor toroidal radii) the equivalent surface tension created by line typing plus finite Larmor radius effects will suppress the instabilities, as shown in Part II, sub-part L;
(ii) by automatic feedback stabilization.
Hitherto feedback stabilization has been demonstrated convincingly only upon certain open-ended devices (linear theta pinches and quiescent or "Q-machine" magnetic mirror devices) which appear to have static structural stability. This is not surprising, for static structural stability of the plasma configuration is necessary for a property known in systems engineering theory as "identifiability", which in turn is necessary for the rational design of feedback stability augmentation or stabilization systems.
For ordinary linear electro-mechanical systems with a finite number of degrees of freedom, static structural stability is equivalent to the requirement that there be no repeated (or "degenerate") complex frequencies and no vanishing real parts of any complex frequencies, while dynamic asymptotic structural stability is equivalent to the requirement that all complex frequencies must have negative real parts. The simultaneous imposition of these two (logically independent) requirements restricts attention to a class of systems sufficiently "well-behaved" to be amenable to routine systems engineering techniques. For example, a rather general procedure for designing a feedback control system applicable to the present plasma configuration is given in "Magnetic Feedback Stabilization in A Tokomak" by Yu. P. Ladikov and Yu. I. Samoilenko, Zuhrnal Technicheskoi Fiziki, Vol. 42, No. 10, pp. 2062-2073, Oct. 1972 (English translation in Soviet Physics--Technical Physics, Vol. 17, No. 10, pp. 1644-1650, April 1973). A similar consideration for nonlinear systems with infinitely many degrees of freedom is implied in the principles underlying the present invention.
The method of computing the external currents required by the principles of this invention is the standard numerical problem in the art of solving elliptic partial differential equations; see, for example, the section on "Toroidal Equilibra With Scalar Pressure" in the chapter on Computational Problems by J. Kileen, pp. 231, 232 of Physics of Hot Plasmas, edited by Rye and Taylor, Oliver and Boyd, 1970. In fact, inserting equations (7) into the vacuum magnetic field equations
.gradient..multidot.B=0 , .gradient..times.B=0 (12)
one finds the well-known Stokes' equation: ##EQU13## which is to be solved outside of .GAMMA. with the Cauchy data (8) and (11) in the manner explained above and illustrated in FIG. 13. The specific numerical examples shown in FIGS. 10, 11 and 12 were computed by applying a standard five-point difference equation approximation to equations (8), (11) and (13).
A. EQUILIBRIUM BOUNDARY CONDITION
The equations for equilibrium in magneto-fluid dynamics are
.gradient..multidot.B=0, (1.1)
.gradient..times.B=j, (1.2)
j.times.B =.gradient.p. (1.3)
Let V denote a volume in the shape of a small "pill box" with one face of area .DELTA.S inside and one parallel to the normal n to the surface .SIGMA. at a point inside the volume. The element of area of the outer face is then approximately dS.congruent.n.DELTA.S while that of the inner face is dS.congruent.-n.DELTA.S. The following calculations become accurate in the limit as the volume V is taken to decrease to zero. By the Divergence Theorem of Gauss, ##EQU14## where (1.1) is used. The surface integral is approximately (after canceling .DELTA.S)
[n.multidot.B].sub.o -[n.multidot.B].sub.i =0 (1.5)
where the subscripts o and i denote outside and inside of .SIGMA., respectively. Thus the normal component of B must be continuous across .SIGMA..
Next, use (1.1) again to obtain the identity
.gradient..multidot.BB=B(.gradient..multidot.B)+(B.multidot..gradient.)B=(B.multidot..gradient.)B (1.6)
pertaining to the dyad BB. Now combine (1.2) and (1.3) to find
.gradient.p=(.gradient..times.B).times.B=(B.multidot..gradient.)B-1/2 .gradient.B.multidot.B=.gradient..multidot.BB-1/2.gradient.B.sup.2
or
.gradient.(p+1/2B.sup.2)-.gradient..multidot.BB=0. (1.7)
Thus, as before
.intg..intg..intg.{.gradient.(p+1/2B.sup.2)-.gradient..multidot.BB}dV=0=.intg..intg.{(p+1/2B.sup.2)n-(n.multidot.B)B}dS=.congruent.({(p+1/2B.sup.2).sub.o -(p+1/2B.sup.2).sub.i }n-{[(n.multidot.B)B].sub.o -[(n.multidot.B)B].sub.1 }).DELTA.S (1.8)
Take the dot product of (1.8) with n and divide by .DELTA.S to obtain
(p+1/2B.sup.2).sub.o -(p+1/2B.sup.2).sub.i -{(n.multidot.B).sub.o.sup.2 -(n.multidot.B).sub.i.sup.2 }=0. (1.9)
Now apply (1.5) in (1.9) to obtain, finally
(p+1/2B.sup.2).sub.o -(p+1/2B.sup.2).sub.i =0, (1.10)
thus proving that the total pressure p.sub.T must be continuous across .SIGMA., where
p.sub.T =p+1/2B.sup.2. (1.11)
B. ANALYTICITY OF THE BOUNDARY CONDITION
In order to establish the physical realizability of the magnetic field configuration of FIGS. 3a and 3b, it is necessary to show that equations (1.1)-(1.3) have a regular (real anayltic) solution satisfying the boundary conditions (1.5) and (1.10) in some external neighborhood of .SIGMA., i.e., that there need be no external conductors tangent to or intersecting .SIGMA. at any point. Indeed, the vacuum magnetic field cannot exist where B has a singularity in the sense of analytic function theory); hence the loci of singularities of B is necessarily a subset of the loci of required external conductors. The nearest external conductor must therefore be at least as close to .SIGMA. as the nearest singularity of B. Hence it is required to show that the nearest singularity of B lies at a finite (non-infinitesimal) distance from the plasma, i.e. that no singularities of B lie on the surface .SIGMA..
This can be accomplished by quoting the Cauchy-Kovalevski theorem. From the physical point of view, since .SIGMA. is freely disposable by the reactor designer, one may without loss of generality assume that .SIGMA. is a real-analytic two-manifold. It will be shown that B is determined by an analytic partial differential equation, so that it remains only to show that the boundary conditions determining B are also analytic in order for the physical realizability of the required magnetic configuration to be established.
This calculation is greatly facilitated by (though it does not necessitate) the use of the assumption of axial symmetry and the subsequent introduction of appropriate cylindrical coordinates (r,.theta.,z).
Under the requirement of axial symmetry (.differential./.differential..theta..ident.0), equation (1.1) becomes ##EQU15## which can be satisfied by choosing ##EQU16## where .PSI. is a function of r and z, as yet arbitrary, called the flux function. Equations (1.2) and (1.3) may be combined to eliminate J, yielding
(.gradient..times.B).times.B=.gradient.p. (2.3)
The .theta. component of equation (2.3) is the Jacobian identity ##EQU17## which implies that ##EQU18## where f is an arbitrary function of .PSI.. Similarly, the other two components of equation (2.3) now give ##EQU19## which implies that
p=g(.PSI.) (2.7)
where g is an arbitrary fucntion of .PSI., and ##EQU20## where ' denotes differentiation with respect to .PSI.. A proof that the same equation holds also for a class of non-zero steady fluid velocities is given in sub-part H.
These equations hold both inside and outside the plasma. For the case of low or medium power production per unit volume (which implies .beta.=1), equation (2.8) must be solved with suitable functions f and g. A choice frequently used in this art is to take f and g to be linear functions, with constant co-efficients, of .psi., in which case equation (2.8) becomes linear in .psi.. The resultant equation is separable and has a well known complete system of eigenfunctions in terms of which every possible solution may be expanded. A more complicated but more accurate procedure would be to choose f and g so that equation (2.8) is compatible with Spitzer's generalization of Ohm's Law (see equation (10.8) below). For ease of illustration, we shall confine further attention principally to the .beta.=1 case, although the principles of the present invention are applicable to arbitrary values of .beta. (0<.beta.<1) in the manner just indicated.
For the case of maximum power production per unit volume (which implies .beta.=1) the equations for the interior of the plasma simply reduce to the equation.
p=constant (2.9)
Thus, on .SIGMA., .psi.=constant. We may, without loss of generality, choose this constant to be zero.
Outside .SIGMA., since p and j both vanish, we see that f'(.PSI.).ident.0 and g(.PSI.).ident.0. Thus, outside .SIGMA. equation (2.8) reduces to Stokes' Equation, ##EQU21## while equation (2.5) becomes ##EQU22## where c is some constant.
Note that .gradient..PSI. is normal to .SIGMA. and that B.multidot..gradient..PSI.=0 as required by equation (1.5) when B.sub.i =0.
We illustrate this situation in FIGS. 4a and 4b. In FIG. 4a, 14 is the topological torus representing the plasma volume; the surface of this is .SIGMA.. The poloidal cross-section 6-7 of this torus is shown in FIG. 4b where 141 represents the cross-section of the plasma volume; .GAMMA. is the poloidal cross-section of .SIGMA. and is the curve
.GAMMA.: .PSI.(r, z)=0. (2.12)
We have illustrated a curve which is symmetircal relative to the r-axis. This is done for convenience, and a magnetic confinement device may be constructed in this way, but it is not an essential requirement.
It is easy to verify that, on .GAMMA., equation (5) of Part 2A above may now be written as ##EQU23## or, equivalently, on .GAMMA.
.vertline..gradient..PSI..vertline.=2.mu.pr.sup.2 -c.sup.2. (2.13)
Clearly, for small enough r, the right-hand side of equation 2.l3 is <0, whereas the left-hand side is always .gtoreq.0. Therefore, there is a minimum r, r.sub.c, at which equation (2.13) can be satisfied, i.e., at which a physically real solution exists. Thus, in FIGS. 4a and 4b, the plasma volume 141 must always lie to the right of the vertical line r=r.sub.c. The critical radius r.sub.c is given by ##EQU24## and we have ##EQU25##
If .GAMMA. lies wholly to the right of r.sub.c, then a particular value of the .+-. sign choice in (2.15) must be made, for .differential..PSI./.differential.n cannot then vanish anywhere on .GAMMA.. But this implies that the vector filed on .SIGMA. necessarily has a non-zero rotational transform angle, which violates requirement (a) of the structural stability criteria of the present method.
It will be seen in sub-parts E and F that .GAMMA. cannot intersect the line r=r.sub.c in an odd or an infinite number of times without violating the structural stability criteria. The only remaining possibility is that .GAMMA. intersect the line r=r.sub.c twice, as in FIG. 4b, or an even number N of times, as illustrated in FIG. 5 for the case N=4.
The problem of physical realizeability now becomes the problem of solving Stokes' Equation
.DELTA..PSI.=0 (2.16)
in an exterior neighborhood of .GAMMA. together with the boundary conditions ##EQU26## on .GAMMA.. Here ##EQU27## In order to define this problem more precisely, the following assumptions will be made.
(i) The boundary .GAMMA. is defined parametrically by r=r(t), z=z(t), 0.ltoreq.t.ltoreq.2.pi., where r(t) and z(t) are real analytic functions of t (i.e. are equal to their Taylor series for each t).
(ii) for all t, 0.ltoreq.t <2.pi., {r}.sup.2 +{z}.sup.2 >0.
(iii) The equation r=r.sub.c holds only at points where the boundary .GAMMA. is nearest to the z-axis, i.e., r(t) is a minimum for t=t.sub.o where r(t.sub.o)=r.sub.c.
(iv) The curvature of the boundary .GAMMA. is not zero when r=r.sub.c.
(v) The indicator .epsilon. changes from +.sqroot.2.mu.p to -.sqroot.2.mu.p or vice versa above and below a point where r=r.sub.c.
Under the preceding assumptions, the function .epsilon. .sqroot.r.sup.2 -r.sub.c.sup.2 is real analytic at all points of the boundary .GAMMA., as will now be shown.
The function .epsilon..sqroot.r.sup.2 -r.sub.c.sup.2 is analytic at points where r>r.sub.c because .epsilon. is a constant in a neighborhood of such a point and the square root of positive real analytic functions is real analytic.
For a point (r.sub.c,z.sub.o) where assumptions (i) through (v) are satisfied, z(t) has a real analytic inverse t(z) in a neighborhood of z.sub.o by the inverse function theorem, where z.sub.o =z(t.sub.o), r.sub.c =r(t.sub.o); indeed, r(t.sub.o)=0 by (iii) so that z(t.sub.o).noteq.0 by (ii).
Thus r=r(t(z))=r(z) is an analytic function of z in a neighborhood of z.sub.o. By (iii), r'(z.sub.o)=0 and, by (iv), r.increment.(z.sub.o)>0. Hence ##EQU28## Both .epsilon. and (z-z.sub.o) change sign at z.sub.o so either .epsilon.(z-z.sub.o)=(z-z.sub.o) or .epsilon.(z-z.sub.o)=-(z-z.sub.o), depending on .epsilon., and these are both real analytic. Also the function under the radical sign in (2.20) is real analytic and positive in a neighborhood of z.sub.o so .epsilon..sqroot.r-r.sub.c is real analytic in a neighborhood of (r.sub.c, z.sub.o).
The function .sqroot.r+r.sub.c is real analytic since (r+r.sub.c) >0 by (iii). Thus .epsilon..sqroot.r.sup.2 -r.sub.c.sup.2 =(.sqroot.r+r.sub.c) .ident..epsilon..sqroot.r-r.sub.c } is real analytic at all points of the boundary. Hence the Cauchy-Kovalevski theorem implies that (2.16)-(2.17) has a real analytic solution in some exterior neighborhood of .GAMMA.. Therefore the class of magnetic field configuration shown in FIGS. 3a and 3b, including such generalizations as are suggested in FIG. 5, is physically realizable.
C. GLOBAL FLUX-LINE PORTRAITS BY CONTOUR INTEGRATION
The problem is to solve the equation ##EQU29## in the region between r>0 and outside of an analytic curve .GAMMA., subject to the boundary conditions .PSI.= on .GAMMA. and ##EQU30## on .GAMMA.. .PSI.(r,z) is to be continuous in r and z on .GAMMA.. F(t) is given by .epsilon.A .sqroot.r.sup.2 (t)-r.sub.c.sup.2 and ##EQU31## Here, .mu..sub.o, p, C, A, r.sub.c are positive constant with C.sub.o >C.sub.i. .GAMMA. is a simple, smooth, closed curve given by r=r(t), z=z(t) where r(t) and z(t) are analytic and periodic of period 2.pi.. .GAMMA.lies to the right of r=r.sub.c and touches this line at exactly two different points, P.sub.1 and P.sub.2. As t increases .GAMMA.is traced in a counter-clockwise direction. Proceeding in this direction from P.sub.1 to P.sub.2, .epsilon.=-1. On the remainder of the curve .epsilon.=+1. .differential..PSI./.differential.n is the outward normal to .GAMMA.. Assume r.sup.2 +z.sup.2 >0 ##EQU32## Assume z.sub.1 .noteq.0 at P.sub.1 and P.sub.2 and the curvature at P.sub.1 and P.sub.2 is not zero. Solution:
.PSI.(r,z)=.PHI.(r,z)-.phi.(r,z)
where ##EQU33## .GAMMA.'must include the same singularities of g(w) inside .GAMMA. but none of the singularities of g(w) outside .GAMMA.. a,b is any point inside .GAMMA.'. r,z is inside the circle .xi..sub.o which is the bi-polar coordinate with a,b as coordinate center which reaches the nearest singularity of g(w).
If .GAMMA.'=.GAMMA. or if .xi. is all inside .GAMMA. then .phi.(r,z)=.PHI.(r,z).
Bipolar Coordinates ##EQU34## .phi.(r,z) is defined for all r,z in a circle about the point (a coth.xi..sub.o, b) of radius acsch.xi..sub.o, where the nearest singularity of g(w) outside .GAMMA. is at r.sub.1,z.sub.1 and .xi..sub.o is a solution of the equation (r.sub.1 -a coth.xi..sub.o).sup.2 +(z.sub.1 -b).sup.2 =(a csch .xi..sub.o).sup.2. Since .phi. is also an analytic function of a and b in .GAMMA..sup.1 we can take a,b anywhere in .GAMMA..sup.1 and find the nearest singularity, (r.sub.2, z.sub.2) and compute the new circle of convergence
(r.sub.2 -a coth.xi..sub.o).sup.2 =(z.sub.2 -b).sup.2 =(a csch .xi..sub.0).sup.2
Since .GAMMA.' can include any point in the half-plane r>0 and still not include the singularities of g(w) outside .GAMMA. we need only this continuation to get to any point on the half-plane (except, of course, the singularities). This does not give the whole plane with one circle, but it does give the stream function near any desired point. If the branch cuts are to the left of the singularities, it should be remembered that the circle .xi..sub.o might be tangent to one of these lines rather than through a singularity. If we integrate in the t-plane we must be sure that we stay in the region which is 1--1 with the w-plane.
Application to a Limacon
The only singularity of a limacon outside of .GAMMA. is at t.sub.1 =.pi.-i 1n b/a: ##EQU35##
z=Iw(t)=e.sup.-.sigma. sin.tau.(b+a e.sup.-.sigma. cost.tau.)
P is the principal part of the branch. It depends on the choice of cut. ##EQU36##
Choose the branch >.tau.>2.pi.
Stokes Equation and Legendre Functions ##EQU37## The Cut:
The cut at w.sub.1 can be chosen in at least 3 different ways:
(a) on the real axis to the left of w.sub.1
(b) on the real axis to the right of w.sub.1
(c) below the z axis at r=w.sub.1
(a) and (b) have the advantage of symmetry with respect to the r axis. (b) has the advantage that the two cuts overlap. To the left of w.sub.1 are no singularities. (a) and (c) have the advantage that the cut doesn't go through the Limacon. In (b) the top and bottom parts of the Limacon are in separate regions of the t-plane. (b) and (c) have the advantage that there is no singularity to the left of w.sub.1. Since we are interested in values on the real axis, (c) seems to be the best choice, despite the lack of symmetry. The 3 cuts are indicated in FIGS. 10-13.
Angles:
In (a) .theta. is chosen so that -.pi.<.theta.<.pi.
In (b) .theta. is chosen so that 0 <.theta.<2.pi.
In (c) .theta. is chosen so that (-.pi.)/2 <.theta.<3.pi./2
The coordinates for case (c) are: ##EQU38##
There are 4 curves in the t plane corresponding to z=0.
(1) .tau.=.pi. .sigma.>-ln b/a
(2) .pi./2<.tau.<.pi. .sigma.=-ln (-b/a COS.tau.)
(3) .tau.=0 -.infin.<.sigma.<.infin.
(4) .tau.=2.pi. -.infin.<.sigma.<.infin.
Jorgens give a real integral ##EQU39## This is equivalent to a line integral on .GAMMA. ##EQU40## where the integral goes from H to F, in a counter-clockwise direction. See FIG. 14.
This is not a closed path unless it avoids the cut as indicated in FIG. 15.
When we distort .GAMMA., to a new curve .GAMMA..sup.1, the beginning point H and the ending point F must not change, if the a.sub.n is to be invariant. See FIG. 16.
The second example would be okay if we included integrals along the cut from H.sup.1 to H and from F to F.sup.1. In the t-plane this is shown in FIG. 17.
If we wanted .GAMMA..sup.1 to cover the whole plane to the right of z=0 and outside .GAMMA. we would integrate from H to F as shown in FIG. 18. ##EQU41## Any contour inside of this contour could be taken if the numerical work is simpler. However, the path .GAMMA..sup.1 should include the segment from A to l if we are inerested in finding flux lines to the left of W.sub.1. In locating a to the right of W.sub..infin. the nearest singularity is the branch cut. The integral (16) can be distorted in a similar way if one avoids the point r,z.
After .PSI. has been computed at every point, the curves for the fluxlines .PSI.(r,z)=constant can be computed from ##EQU42## r(s.sub.0)=r.sub.0, z(s.sub.0)=z.sub.0 This gives r=r(s), z=z(s) as the parametric equations of the flux lines.
D. NUMERICAL COMPUTATION OF SOLUTIONS OF STOKES' EQUATION WITH PRESSURE BALANCE BOUNDARY CONDITIONS
INTRODUCTION
We describe the results of applying difference methods to the numerical solution of an elliptic partial differential equation (Stokes') in the region exterior to a non-convex boundary which boundary has Cauchy-type conditions involving both values of the solution and of the normal derivative. Five point difference equations (Laplace type) corresponding to Stokes' equation are constructed in the usual manner for a conformal but noncartesian grid. Given the value of the function on the boundary and a second order approximation of the function via the normal derivative in the first exterior difference layer off of the boundary, the function can be continued into the exterior by employing the five point difference scheme. Although this continuation cannot be carried out indefinitely, the solutions we seek are essentially local ones; that is, we usually need only to go beyond one singularity outside the plasma boundary. FIG. 19 is a sketch summarizing the numerical results and establishes our notation. The subsequent FIGS. 21-23 are scale reproductions of actual calculated plots and upon fixing a coordinate system may be regarded as numerical data suitable for the construction of specific machines. (The dotted curves are the expected extrapolations beyond the range of the difference procedure.)
A DIFFERENCE METHOD
The purpose of this section is two-fold. One is to locate a reasonable source of non-convex regions of the type suggested by FIG. 19. We should thus be able (without seriously altering the numerical technique for solving Stokes' equation for such a boundary) to control the aspect ratio, dimple, etc., of the region. The second is to develop the difference equations germane to such a family.
the boundary value problem exposed above consisted of a simple closed analytic curve, .GAMMA., lying in the right half plane. The interior region associated with this curve is assumed to abe a non-convex region as in FIG. 19 (or its mirror-image in the y-axis) and touches the critical line at exactly two points. We assume for convenience that the region is symmetric with respect to the horizontal axis. We suppose that the solution, .PSI., of Stokes' equation is zero on the curve .GAMMA. and its (outwardly directed) normal derivative equal to a prescribed function derived from the pressure balance of the plasma and the magnetic field at the plasma surface. (We view .GAMMA. as a cross-section of this plasma surface.) The problem is to solve the system
(a) .DELTA..sub.* .PSI.=0 (exterior to .GAMMA.)
(b) .differential..PSI./.differential.n=f (equality on the curve .GAMMA.)
(c) .PSI.=0 (equality on the curve .GAMMA.)
where .PSI.=.PSI.(x,y), f=.epsilon..sqroot.x.sup.2 -x.sub.c.sup.2,
the values of .epsilon. are indicated in FIG. 19, n is the (outwardly-directed) normal, .differential./.differential.n is the normal derivative and .DELTA. is the (elliptic) Stokes' operator ##EQU43## We refer to this as a Cauchy-Stokes system.
Since .GAMMA. is a simple closed curve, one expects that a reasonable source of such curves can be had by conformal mapping of a line segment. Furthermore, a line segment has a natural grid associated with it, hence the conformal map will carry this grid into a natural grid around the curve .GAMMA.. We will compose several conformal maps as follows. Refer to the z-plane with coordinates x,y,z=x+iy, and the W-plane u,v,w=u+iv. Take the line segment, -.pi.<x.ltoreq..pi., y=0 as in FIG. 20a, map to the W-plane in FIG. 20b, via the sin z map. Note that the ellipse is the image of the line segment y=y.sub.o, -.pi.<x.ltoreq..pi.. Then rotate the ellipse 90.degree. and translate to the right (this translation relates to the aspect ratio of the torus). Then the right half of the plane in FIG. 20d is the result of squaring a hyperbolic region of the right half plane in FIG. 20c. Of course, this half plane in FIG. 20d containing the nonconvex region is the domain of interest for our solutions of Stokes' equation. A specific example of this procedure is the conformal map ##EQU44## where the z refers to the plane in FIG. 20a and the w to the plane in FIG. 20d. u.sub.0 is the distance of the center of the ellipse to the origin in the plane of FIG. 20c; u.sub.1 is the distance from the singularity to the origin after the region of FIG. 20d is translated for a particular aspect ratio.
We now proceed to the development of Stokes' analogues of the Laplace five point formula. We will solve the difference equation in the z-plane which yields, under the conformal map, the solution of Stokes' equation in the w-plane. Naturally the two boundary conditions to be imposed in the z-plane will differ from those of the original problem. It is convenient to let .PSI.=x.sup.1/2 Q, which gives rise to the Cauchy-Stokes system in the z-plane: ##EQU45##
Converting this system by the usual procedure into a difference system yields the five point equation (corresponding to (a) above) ##EQU46## where h and k are independent increments in the y and x directions respectively. The (c) boundary condition simply states that Q is zero for the boundary points, viz., the first layer of the difference scheme. The second layer obtains by applying condition (b) and recognizing that Q has (locally) a Taylor expansion along the direction of the normal to .GAMMA.. (If s is direction along the normal then ##EQU47## represents three terms of that expansion.) Thus conditions (b) and (c) give the first two layers and the five point difference formula is operative to get the third layer.
THE NUMERICAL PLOTS
This section consists of the pictures of a number of numerical experimentsusing the techniques of the previous section. These involve variations of distance of the plasma from the origin, height and width of the plasma, dimple size, as well as variations of the ratio of h and k and mesh size for observing the good stability (numerically) of the difference scheme. The computer drawn pictures present three flux lines, viz., .PSI.=0 and .PSI.=two other constants, one line around the singularity and the other line around the plasma, singularity and .PSI.=0 line. (We define a flux line to be a .PSI.=constant curve.) Note that the pictures in this section are of .PSI. and not of Q curves. u.sub.1 and u.sub.0 are defined in the previous section; .alpha.=h/k and t=b/a where b is the minor and a is the major axis of the ellipse in FIG. 20c. FIG. 21 includes the conformal grid associated with the difference scheme. (For consistency with FIGS. 1 and 2, FIGS. 21, 22 and 23 have been mirror-reflected in the y-axis, after being originally computed precisely as indicated above.)
The parametric form of the equation for the curve, .GAMMA., describing the poloidal plasma cross-section for the parameters of FIG. 21 is
.omega.=25.1503125+0.225 cos .chi.+0.1503125 cos 2.chi.+ i[0.5 sin .chi.+0.1125 sin 2.chi.]
where the parameter .chi. varies between -.pi. and .pi.. Any units may be used for .omega.. In FIG. 21, 0.45 .omega.-units is 2 inches. The toroidal axis is 9.454 feet to the right of the outer extremity of the plasma. The line of symmetry is the radius axis.
FIG. 22 shows the somewhat enlarged plasma surface, .GAMMA., and flux surfaces for the same parameters as for FIG. 21 but with the grid removed so as to exhibit more clearly the nature of the flux surfaces.
FIG. 23 is similar to FIG. 22, but with a different set of parameters.
These numerical plots exhibit clearly the physical realizability of the structurally stable magnetic field configuration of this invention. These magnetic fields may be practically realized in a device by locating current carrying wires along a given flux surface and everywhere that point and line singularities of these flux surfaces occur; or locating current double layers (see under Part I above) along surfaces which are not flux surfaces in accordance with standard procedures in this art (see, e.g., Stratton, Electromagnetic Theory, McGraw-Hill, 1941). FIG. 2B shows an example of this latter method. Since the magnetic configuration described corresponds to a non-zero net flux through the surface between the plasma and the toroidal axis, an equal amount of bias flux must be initially present passing somewhere through this same surface, preferably near the toroidal axis. This bias flux must be made to vanish as the confining magnetic configuration is established.
Alternatively, the net flux through the surface between the plasma and the toroidal axis may be kept at zero by turning on a biasflux near the axis, equal in magnitude but opposite in direction, to the net flux through this surface due to the main field as the main field is turned on. E. NUMBER OF LIMIT CYCLES MUST BE EVEN
According to criterion (d) of structural stability specified in Part 2A above, the unit-length vector field B/B on .pi. can have only a finite number N of closed orbits (i.e. periodic solutions of dr/ds=B(r)/B(r)); and each such orbit must be simple in the sense that equation (6) must hold. But equation (6) implies that (see [5.1]) each such orbit must be asymptotically stable either for the given sense of the arc length parameter s, or after s is replaced by (-s). This clearly rules out the possibility of a semi-stable limit cycle. But if the number N were odd, there would have to be at least one limit cycle with trajectories approaching it on one side and leaving it on the other side, i.e., a semi-stable limit cycle. This contradiction proves that N must be an even number.
REFERENCE
[5.1] S. Lefschetz, Differential Equations: Geometric Theory, Inter-Science, 1963.
F. ROTATIONAL TRANSFORM ANGLE MUST BE ZERO
Here we employ the Poincare-Denjoy theory of differential equations on a torus. See, for example, Kneser [6.1], Coddington-Levinson [6.2], and hartman [6.3]. In the standard mathematical theory the rotation number lies between zero and unity because the arc length of the appropriate circle of minor cross-section is normalized to unity. In the present art the rotational transform angle is the same as the standard mathematical rotation number multiplied by 2.pi..
We now restrict attention to the magnetic field lines on a toroidal magnetic shell which is the surface .SIGMA. of a unity .beta. plasma. The normal component of B must be continuous across .SIGMA., and since there is no magnetic field inside of .SIGMA., the normal component must also vanish outside of .SIGMA., i.e., 93 must be tangent to .SIGMA.. Also the combination of Maxwell's Equations of electromagnetism and Bernoulli's principle for fluids implies that the total pressure p.sub.T must be continuous across .SIGMA., i.e. B=.sqroot.2.mu.p>0. Hence the vector .SIGMA. not only cannot vanish on .SIGMA. but it must have a constant length B>0 everywhere on .SIGMA..
By hypothesis the toroidal component B.sub.0 of B is always greater than the poloidal component .sqroot.B.sub.r.sup.2 +B.sub.z.sup.2 (assuming axial symmetry). Hence the vector B always cuts the topological circle .GAMMA. at a non-zero angle (i.e. it cannot be tangent to .GAMMA.). Therefore if r.sub.o is an arbitrary point on .GAMMA., the solution r=r(s) of dr/ds=B(r)/B(r) which satisfies r(0)=r.sub.o must have a smallest value S.sub.1 >0 for which r.sub.1 =r(S.sub.1) lies on .GAMMA.; likewise, there will be a smallest S.sub.k+1 >S.sub.k for which r.sub.k+1 =r(S.sub.k+1) again lies on .GAMMA., for k=0,1,2,3, . . . Let .delta..sub.k denote the arc length on .GAMMA. between r.sub.k and r.sub.k+1 divided by the total arc length of .GAMMA.. The rotational transform angle .alpha. is defined by ##EQU48## It is proved in the works cited that .alpha. is independent of the particular choice of r.sub.o on .GAMMA.. The conventional mathematical rotation number is now (.alpha./2.pi.).
For the present purpose it is reasonable to say that, by definition, an equilibrium configuration (.SIGMA..sub.o,B.sub.o) lacks static structural stability if there is an arbitrarily close configuration (.SIGMA..sub.1,B.sub.1) which lacks structural stability in the mathematical sense of Peixoto and Pliss and which is physically accessible from (.SIGMA.,B), i.e. which can be attained by a continuous deformation (.SIGMA..sub..sigma., B.sub..sigma.), 0.ltoreq..sigma..ltoreq.1, in which each configuration (.SIGMA..sub..sigma.,B.sub.94 ) is a magneto-static equilibrium which is physically realizeable for some plasma pressure p.sub..sigma. >0.
The necessary and sufficient criteria for structural stability stated by Pliss are equivalent to those of Peixoto but are formulated somewhat differently; they are:
(a) the vector field must have a rational rotation number;
(b) all the periodic field lines must have non-zero characteristic exponents.
The rotational transform angle has been used in the present art by Spitzer, (U.S. Pat. No. 3,002,912), by Kulsrud et al (U.S. Pat. No. 3,012,955) and by Koenig (U.S. Pat. No. 3,088,894), wherein magneto-static equilibrium of the plasma is achieved by means of a non-zero rotational transform. It will be shown below that if one reduces the rotational transform in a Stellarator, Tokomak, Screw Pinch, Scyllac (or high-.beta. Stellarator), or similr device toward zero then of necessity one must bring the plasma into contact with an external conductor, obviating the device's utility.
In the simplest case, .GAMMA. may be taken to be the boundary of a convex body, e.g. an ellipse or a circle, and there turns out to be a simple dichotomy between the cases where the rotation number is a rational number (terminating or repeating decimal fraction) and the cases where it is an irrational number (infinite non-repeating decimal fraction): when (.alpha./2.pi.) is rational all of the field lines are closed, and when it is irrational all of them are "ergodic" (i.e. wind densely on the torus .SIGMA. without ever exactly closing upon themselves).
In the present art it is seemingly universally believed and stated, erroneously, that the preceding dichotomy is true in general. For example, the 1965 survey article by Greene and Johnson ([6.4]) repeats this fallacious dichotomy on page 201 (lines 8-9 and 20-21) and then uses it as a basis for the remainder of the 43-page survey. However, there is in actuality a third possibility, treated definitively in 1924 by Kneser [2.1], in which all but a finite number of the field lines are "asymptotic" (to one of finite number of closed limit cycles). It is this third case upon which attention is focused by Pliss's structural stability criteria.
Now consider the case of axial symmetry. If .GAMMA. lies wholly to the right of the vertical line r=r.sub.c, then .differential..PSI./.differential.n.noteq.0 at all points on .GAMMA., and a particular choice of the .+-. sign in equation (2.15) must be made; for example, choose the + sign. Now .GAMMA. may be chosen to be an analytic curve without essential loss of generality, and since .PSI.=0 and .differential..PSI./.differential.n=C.sqroot.r.sup.2 r.sub.c.sup.-2 -1 are (for r>r.sub.c) analytic boundary conditions on .GAMMA. the Cauchy-Kovalevski Theorem implies that .PSI. is a real analytic function of (r,z) in a neighborhood of .GAMMA. (increasing outward from .GAMMA.). Thus in a neighborhood of .GAMMA. the curves .PSI.=constant form a set of nested circles which by rotation about the axis of symmetry generate a family of nested tori about .SIGMA.. Now (.alpha./2.pi.) will be an analytic function of the indexing parameter .PSI.; if not identically constant it will have both rational and irrational values arbitrarily near to .PSI.=0, a case eliminated by the definition of static structural stability. Likewise the case where (.alpha./2.pi.) is an irrational constant can be a priori eliminated from consideration by Pliss's first criterion. The case where (.alpha./2.pi.) is hypothesized to be a rational constant may be eliminated because of the requirements of isolated closed field lines and axial symmetry.
Therefore .GAMMA. must intersect the line r=r.sub.c. No point on .GAMMA. may lie to the left of this line (i.e., in the region r>r.sub.c) because then the square root in equation (2.15) would have to be imaginary. Thus .GAMMA. must touch the line r=r.sub.c but not cross it.
At each point where .GAMMA. touches the line r=r.sub.c, i.e., where .differential..PSI./.differential.n=0, one must have B.sub.r =B.sub.z =0. Hence the field line through such a point is purely azimuthal and, by axial symmetry, continuous around the torus, and forms a closed orbit.
Hence .GAMMA. can touch the line r=r.sub.c only at a finite number of points, by Peixoto's criterion (d), and all other field lines or "trajectories" must be asymptotic to these closed orbits by Peixoto's criterion (b); hence these closed orbits are limit cycles. As above in sub-part E, there can only be an even number of them.
We shall next show that in this case the rotational transform must be identically zero. Pick a point on .GAMMA. whose arc-length distance from its limit-cycle's intersection with .GAMMA., divided by the arc length of .GAMMA., is less than a pre-specified arbitrary constant .epsilon.. Now the trajectory through this point spirals endlessly nearer to the limit-cycle without crossing it, in such a manner that each .delta..sub. defined above satisfies .delta..sub.k =.epsilon.. Hence ##EQU49## But .epsilon. was arbitrary and (.alpha./2.pi.) <.epsilon. is independent of the original point chosen on .GAMMA.. Hence (.alpha./2.pi.)=0.
Thus the requirement that .GAMMA. must touch the line r=r.sub.c exactly an even number of times, as illustrated in FIG. 5, implies that the rotational transform .alpha. of the field on .SIGMA. must vanish.
The rotational transform angle .alpha. may be calculated explicitly for the case where .GAMMA. is a circle, as in Jorgens [6.5] (see also Lust and Schluter [6.6], Biermann et al [6.7], and Kippenhahn [6.8]). Letting .alpha. tend to zero implies that .GAMMA. moves closer to r=r.sub.c, finally touching the line at the point 18 in FIG. 9 when .alpha.=0. The results of Jorgens just cited imply that point 18 is an elliptic singularity and that there must be an external conductor through point 18. Hence the usual case where .GAMMA. is the boundary of a convex body, i.e. where the poloidal cross-section of the plasma is convex, cannot be physically realizeable if the rotational transform is required to vanish.
If, however, the square root in equation (2.15) changes sign at each of the even number of points where .GAMMA. touches the line r=r.sub.c then, by the preceding, the rotational transform vanishes identically, while in sub-part B, the corresponding configuration is physically realizable.
In this case the plasma surface cannot leave the line r=r.sub.c. There will be an even number of oppositely-directed toroidal surface currents in strips and therefore (by continuity) there must be an even number of intervening zero-current lines or strips (in between the current strips). But where the toroidal current is zero the poloidal field is also zero and so pressure balance requires that such points touch the line r=r.sub.c. Therefore the hypothesis of a continuous infinitesimal movement of .GAMMA. into the open set r>r.sub.c contradicts the physical existence of the postulated types of sheet currents. One might consider a configuration topologically equivalent to that under consideration except that one or more of the original limit-cycles would be replaced by a continuum of closed orbits; however for infinitesimal perturbations this would violate the Peixoto-Pliss criteria.
REFERENCES
[6.1] H. Kneser, Mathematische Annalen, vol. 91 (1924), pp. 135-154 [6.2] E. A. Coddington and N. Levinson, Theory of Ordinary Differential Equations, McGraw Hill, 1955.
[6.3] P. Hartman, Ordinary Differential Equations, Wiley, 1964.
[6.4] K. Brueckner, editor, Advances in Theoretical Physics, vol. 1 (1965).
[6.5] K. Jorgens, Zeitschrift fur Naturforschung, vol. 13a (1958), pp. 493-498. (English translation available from McDonnell-Douglas Astronautics Company, Santa Monica, California, as paper MDC-G1140.)
[6.6] R. Lust and A. Schluter, Z. Astrophys., vol. 38 (1955), p. 190, and Zeitschrift fur Naturforschung, vol. 12a (1957), p. 850.
[6.7] L. Biermann, K. Hain, K. Jorgens, and R. Lur Naturforschung, vol. 12a (1957), p. 826.
[6.8] R. Kippenhahn, Zeitschrift fur Naturforschung, vol. 13a (1958,), p. 260.
[6.9] J. M. Greene and J. L. Johnson, Physics of Fluids, vol. 4 (1961), pp. 875-890.
G. LIMIT CYCLES ARE IN FACT SIMPLE
It was proved above in sub-parts B, E and F that for static structural stability the plasma boundary's cross section must touch the critical line r=r.sub.c an even number of times and that if the square root in equation (2.15) changes sign at each point of contact the corresponding plasma equilibrium is physically realizable. It remains only to prove that each such contact point generates a "simple" limit cycle in order to fulfill all of the criteria of Peixoto and Pliss.
Here the following topics relative to a vector field on a surface will be treated:
(1) "Area-preserving on a surface" will be defined.
(2) It will be shown that "area-preserving on a surface" is equivalent to the vanishing of the two-dimensional divergence (which quantity will also be defined).
(3) It will be shown that the magnetic field in the method described in Part 2A cannot be area-preserving on the plasma surface. (Its two-dimensional divergence cannot vanish everywhere on the surface).
(4) It will also be shown by means of local solution of the equations that the two-dimensional divergence of the magnetic field cannot even vanish on the limit cycles. It is shown finally that this proves the simplicity criterion of Peixoto and Pliss.
We turn now to a more detailed consideration of the preceding four points:
(1) Consider a surface .SIGMA. with a tangent vector field B. (In the application to the method described in Part 2A, .SIGMA. will be the plasma surface and B will be the magnetic field). Take a closed oriented curve C lying in .SIGMA. and let the vector element of length be dx. Let n be the unit normal vector on .SIGMA. (into the paper). Then n.dx=0 and n.times.dx is a vector of length .vertline.dx.vertline., tangent to .SIGMA., and perpendicular to and pointing out from C. See FIG. 24.
We say that B is area-preserving on .SIGMA. if ##EQU50## for every C lying in .SIGMA.. This says that there are just as many B-lines going out of C as coming in. This is the exact analogue of ##EQU51## for a closed surface .sigma., which is the condition for B to be volume-preserving and which is equivalent to .gradient..B=0 (the three-dimensional divergence vanishes).
(2) We now use the techniques of differential forms to define the two-dimensional divergence of B, .gradient..sup.(2).B, and to prove the equivalence stated above in (2).
We consider two-unit vectors e.sub.1 and e.sub.2 which span the surface .SIGMA.. The vectors e.sub.1,e.sub.2, and n may be taken as a positively-oriented orthonormal triple. A displacement dx in .SIGMA. may be written
dx=.omega..sup.1 e.sub.1+ .omega..sup.2 e.sub.2
where the .omega..sup.i =e.sub.i .multidot.dr are 1-forms. (For background, see H. Flanders, Differential Forms, Academic Press, 1963, pp. 40-44.) The star operator * is defined by requiring it to be linear and by the equations
*.omega..sup.1 =.omega..sup.2
*.omega..sup.2 =-.omega..sup.1
*1=.omega..sup.1 .omega..sup.2
*(.omega..sup.1 107 .sup.2)=1
If a vector B tangent to .SIGMA. is written
B=B.sub.1 e.sub.1 +B.sub.2 e.sub.2
we define a 1-form by
.omega.=B.sub.1 .omega..sup.1 +B.sub.2 .omega..sup.2.
The two-dimensional divergence of B is now defined to be
.gradient..sup.(2) .multidot.B=*d(*.omega.).
(This coincides with the usual expression in Cartesian coordinates.) We note that ##EQU52## Also, d(*.omega.)=*(.gradient..sup.(2) .multidot.B)=.gradient..sup.(2) .multidot.B.omega..sup.1 .omega..sup.2. We integrate the latter equation over the interior S of a closed curve C and use Stokes' theorem: ##EQU53## If now B is area-preserving on .SIGMA., then the last integral vanishes for all C and so .gradient..sup.(2) .multidot.B=0. The converse holds, and so the two conditions are equivalent.
(3) We calculate .gradient..sup.(2) .multidot.B on the surface .PSI.=constant (which is .SIGMA.) for the method described in Part 2A. We note first that ##EQU54## where subscripts denote partial derivatives and r and k are unit vectors in the r and z directions in cylindrical polar coordinates. We take e.sub.1 =.phi.. Then ##EQU55## We note that on .SIGMA., d.PSI.=0, or
.PSI..sub.r dr+.PSI..sub.z d.sub.z =0.
Now
x=rr+zk,
whence
dx=drr+rd.phi..phi.+dzk,
.omega..sup.1 =e.sub.1 .multidot.dx=rd.phi.,
and ##EQU56## Next ##EQU57## so that ##EQU58## Now ##EQU59## Since (r,z)=constant on .SIGMA., we may take r to be a function of z, so that the coefficient of dz is a function of z. The exterior derivative of this term therefore vanishes. Thus, we have ##EQU60## We replace dr, dz, d.sub.100 by their expressions in terms of .omega..sup.1 and .omega..sup.2 and note that .gradient..sup.(2) .multidot.B=*d(*.omega.). We get ##EQU61##
If we wish B to be area-preserving, we set .gradient..sup.(2) .multidot.B=0. This equation can be written in the form
.PSI..sub.z (.PSI..sub.r.sup.2 +.PSI..sub.z.sup.2).sub.r =.PSI..sub.r (.PSI..sub.r.sup.2 +.PSI..sub.z.sup.2).sub.z,
which can be integrated to give
.PSI..sub.r.sup.2 +.PSI..sub.z.sup.2 =F(.PSI.),
where F is an arbitrary function of .PSI.. But on .SIGMA., ##EQU62## so we have ##EQU63## Thus the curve .PSI.=constant must be r=constant, which is impossible for the given design of .SIGMA.. Thus B cannot be area-preserving on .SIGMA., as was to be proved. (4) We construct a solution of equations (2.10), (2.12), (2.15) valid in the neighborhood of a limit-cycle. We denote the position of the limit-cycle as r=r.sub.c, z=z.sub.c, and write .rho.=r-r.sub.c, .zeta.=z-z.sub.c. See FIG. 25. Thus the limit-cycle is at 0 in .rho..zeta. space. Since ##EQU64## at 0, on .SIGMA., we assume an equation for .SIGMA. near 0 as
.rho.=.alpha..zeta..sup.2 +0(.zeta..sup.3)
where .alpha. is a constant. We also assume for .PSI. near 0:
.PSI.=a+b.rho.+h.zeta.+e.rho..sup.2 =f.rho..zeta.+g.zeta..sup.2 +0(3rd order terms).
Equation (2.16), which is now ##EQU65## yields
b=2r.sub.c (g+e).
On .SIGMA., we get
.PSI.=a+b.alpha..zeta..sup.2 +h.zeta.+g.zeta..sup.2 =0(.zeta..sup.3)
.PSI..sub..rho. =b +f.zeta.+0(.zeta..sup.2)
.PSI..sub..zeta. =h+2g.zeta.+0(.zeta..sup.2),
So equation (2.12) gives a=h=0, g=-b.alpha. and equation (2.15) gives b=0 and f.sup.2 r.sub.c =z.alpha.c.sup.2. Thus ##EQU66## (neglecting second order terms). Now .rho. is of order .zeta..sup.2 on .SIGMA., so B.sub.r is small; also, B.sub.100 =c/r.perspectiveto.c/r.sub.c. Thus we get, on .SIGMA. near 0, ##EQU67## which is portrayed in FIG. 26.
The field lines satisfy the equations ##EQU68## where ##EQU69## These integrate to ##EQU70## Finally, note that ##EQU71## at 0, so that ##EQU72## around the limit cycle is ##EQU73## and the limit cycle is simple. Clearly, this applies to any limit cycle on the surface of .SIGMA. when the plasma equilibrium is constructed as above in sub-part B. Therefore all of the criteria of Peixoto and Pliss have been verified and the demonstration that each member of the class of magneto-static plasma equilibria delineated above possesses static structural stability is complete.
H. STATIONARY AZIMUTHAL FLOW
The EMTFD equations at equilibrium with transport coefficients set equal to zero with .eta.=.eta..sub.0 and v.noteq.0 become
.gradient..rho..multidot.v=.rho..gradient..multidot.v=0 (8.1)
-p.gradient..multidot.v=0 (8.2)
(v.multidot..gradient.)v+.gradient.p=(.gradient..times.B).times.B (8.3)
.gradient..times.(v.times.B)=0 (8.4)
Equation (8.2) implies p=0 or .gradient..multidot.v=0. Taking .gradient..multidot.v=0 yields (since p=R.rho..eta..sub.0)
.gradient..rho..multidot.v=0 or .gradient.p.multidot.v=0 (8.1a)
.gradient..multidot.v=0 (8.2a)
(v.multidot..gradient.)v+.gradient.p=(.gradient..times.B).times.B (8.3a)
(B.multidot..gradient.)v=(v.multidot..gradient.)B. (8.4a)
Assuming axial symmetry implies .differential.p/.differential..theta.=0 and with equation (8.2a) implies ##EQU74## .gradient..multidot.B=0 implies ##EQU75## where .PSI. is the magnetic flux function, V is the velocity stream function.
Equation (8.1a) becomes ##EQU76## so p=g(V) if .differential.V/.differential.r or .differential.V/.differential.z.noteq.0.
The first component of (8.4a) becomes ##EQU77## The third component of (8.4a) becomes ##EQU78## These last two imply that ##EQU79## The second component of (8.4a) becomes ##EQU80## The first component of (8.3a) becomes ##EQU81## The third component of (8.3a) becomes ##EQU82## where ##EQU83## The second component of (8.3a) becomes ##EQU84##
In the special case v.sub.r =v.sub.z =0, v.sub..theta. =(const.) r, the above equations imply that V=const. and ##EQU85## where
.DELTA..sub.* .PSI.+f(.PSI.)f'(.PSI.)+r.sup.2 g'(.PSI.)=0.
Hence the analysis in sub-part B is immediately applicable.
I. SINGULARITIES OF LIMACON AND CAP-CYCLIDE CROSS-SECTIONS
Consider a torus of known cross-section .GAMMA.. Let .GAMMA. be given by a simple closed smooth curve,
r=r(.tau.) z=z(.tau.) 0.ltoreq..tau.<2.pi..
Assume r(t),z(t) are analytic functions with two minimum values of r called critical points. Consider the region, R, of the (r,z) plane outside .GAMMA. and with r>0. Let the flux function, .PSI.(x,y), satisfy ##EQU86## and .PSI.=O on .GAMMA. ##EQU87## where c is a positive constant
r.sub.c is the minimum value of r on .GAMMA.
.epsilon. is 1 on .GAMMA. between critical points
.epsilon. is -1 on the remaining part of .GAMMA..
In Jorgen's paper (Zeitschrift fur Naturforschung 13a, 493-498, 1958 "Axially Symmetric Solutions of the Magnetohydrostatic Equations with Surface Currents"), it is shown that the singularities of ##EQU88## include the singularities of the solution of the above problem.
w=r(t)-iz(t)
w=r(t)+iz(t)
t=.tau.+i.sigma.
r(t) and z(t) are the complexification of .GAMMA.. g(w) has singularities when ##EQU89## or r(t)=0,r.sub.c,-r.sub.c, or when w,w are not one-to-one mappings to the t plane. We are interested only in those singularities which are in R.
Case 1: A limacon
W=r.sub.o +(Q+.mu..sub.o).sup.1/2
(See FIG. 27) ##EQU90##
If Q=Re.sup.i.theta. then ##EQU91##
In the w plane ##EQU92##
In the Q plane ##EQU93## In the t plane ##EQU94##
The only singularity of g(w) outside the limacon is when w=0. This is when dw/dQ=0 and is the point A in FIGS. 27A, 27B and 28.
Case 2: Cap Cyclide
The Cap cyclide coordinate system is given by Moon and Spencer in the "Field Theory Handbook". ##EQU95##
The region outside .GAMMA. in the w plane maps into the rectangle indicated. The inside of a circle of radius R in the w plane maps into upper half of the rectangle in the Q plane.
w=r+iz
Q=.mu.-iv
S=sn.mu.=sin .phi. ##EQU96##
Jacobi elliptic functions
S.sub.1 =sn v with k replaced by 1-k.sup.2 =k.sub.1, 0<k.sub.1 <1 ##EQU97## where ##EQU98##
The outside of .GAMMA. in the w plane maps into the inside of a rectangle in the t-plane
t.sub.A =k'/2
t.sub.L =0
t.sub.D -K'/2
t.sub.E =-K'
t.sub.F =-K'-i.mu..sub.o ##EQU99##
t.sub.j =-i.mu..sub.o
t.sub.N =K'-i.mu..sub.o
t.sub.p =K'
Points corresponding to .tau..fwdarw.-.tau. are inversion points in the w plane.
The singularities are as follows:
__________________________________________________________________________Source -u = u.sub.o + .sigma. V t values Comments__________________________________________________________________________.about.w = .infin. 2K - 2uo ##STR1## ##STR2## ##STR3##.w = 0 K K' u.sub.o - it.sub.1 = K + iK K 0 u.sub.o - it.sub.2 = K w.sub.1,w.sub.2 are inside .GAMMA. on circle at ends of curve u = K, w.sub. 2 = -w.sub.1 w.sub.3, w.sub.4 are on circle.about..w = 0 .vertline.K - 2u.sub.o .vertline. 0 u + it.sub.3 = K ##STR4## .vertline.K - 2u.sub.o .vertline. K' u + it.sub.4 = K + iK' w.sub.4 = -w.sub.3r = 0 K - u K' t.sub.5 = K' + iK w.sub.5, w.sub.6 are on circle w.sub.5 = inside if 0 < u.sub.o < K/2 ##STR5##r = .infin. K - u 0 t.sub.6 = iK w.sub.6 = -w.sub.5r = .+-. r.sub.c Solution of w.sub.7, w.sub.8, w.sub.9, w.sub.10 8th degree equation__________________________________________________________________________
Critical points occur at ##EQU100##
If 0<k<3-2 .sqroot.2=0.1716, there will always be a dimple.
If 1>k>3-2.sqroot.2, there will be a dimple only if ##EQU101##
This occurs when z.sub.c =0
If 2/3K<.mu..sub.1, there are 5 singularities outside the plasma ##EQU102##
J. ELECTROMAGNETOTHERMOFLIUD-DYNAMICAL ANALYSIS
The preceding discussion has shown how to synthesize and analyze an idealized magneto-static plasma equilibrium, such as has been regularly employed in the present art for the preliminary design of magnetic confinement devices. In reality when the non-conservative aspects of an ionized gas (including electrical resistivity, thermal conductivity, and hydrodynamic viscosity) are considered, there is no equilibrium unless one assumes a constant influx of mass to compensate for that lost by cross-field diffusion and a steady volume-distributed source of (thermonuclear) energy to compensate for radiative and convective losses.
Moreover there are many phenomena of "micro-instabilities" which may be considered for low-density plasmas.
Recall that the important parameter beta designates the ratio of internal plasma fluid pressure to external confining magnetic field pressure, i.e.
.beta.=p/(B.sub.o.sup.2 /2.mu.) (10.1)
By equation (1.10) it is evident that at the plasma boundary ##EQU103## whence follows the fundamental relation
0<.beta..ltoreq.1. (10.3)
For convenience we shall adopt the following terminology:
______________________________________Low Beta 0 < .beta. < 0.10Medium Beta 0.10 .ltoreq. .beta. < 0.70High Beta 0.70 .ltoreq. .beta. .ltoreq. 1.0______________________________________
The Scylla IV open-ended theta pinch at Los Alamos reported .beta.=0.8 in an experiment of December 1969. The Centaur cusp-ended theta pinch at the U. K. Culham Laboratory has, for hot, dense plasmas, reported values of .beta. in excess of 0.99. Hence it is realistic to analyze high beta plasmas, although nearly all prior work in this art has concerned low beta plasmas (Stellarator, Tokomak, Doublet, etc.).
The preceding results in sub-parts A through I apply chiefly to the vacuum magnetic field, specifically to what in the present art is termed a "magnetic bottle" or "magnetic trap", and are valid regardless of whether the plasma confined in this "bottle" is of low or high beta.
For a given external magnetic field strength, and plasma temperature one may refer to low density or high density plasmas instead of low or high beta, because in this case beta is proportional to the density.
It is customary to examine low density plasmas by mathematical techniques from the kinetic theory of plasmas (see e.g. references [10.1] and [10.2]). For the present purposes, however, it will be adequate to treat low and medium beta plasmas by means of a continuum mechanical model in which the various physical effects (e.g. Hall current) due to finite ion and electron gyro-radii are included by means of a generalization of Spitzer's extended Ohm's Law [10.3]. It will be explained that if the resultant mathematical model of the plasma equilibrium possesses dynamic structural stability, then the neglect of microscopic effects from plasma kinetic theory is justifiable whenever these effects are "sufficiently" small. It will then be shown that the magnitude of these effects decreases to zero as the plasma's beta increases to unity. Hence for high beta the following mathematical model will have assured adequacy for dynamic stability analyses.
An alternative justification for the use of a continuum mechanical model may be obtained from a consideration of the following four expert opinions:
(1) "While high-beta instability theory is not yet sufficiently advanced to answer this question unequivocally, there is reason to believe that micro-instabilities may be weak or possibly even non-existent. If so, this will be another important advantage of high-beta systems." Amasa S. Bishop (Director of Fusion Research, Division of Research, USAEC), lecture of Nov. 12, 1969, before the American Physical Society on "Recent World Developments in Controlled Fusion", published in Nuclear Fusion, vol. 10 (1970), pp. 85-93.
(2) "Although the emphasis in toroidal investigations has long been on microinstability . . . , it has, somewhat belatedly, become clear that these mechanisms should defer in priority to the more basic questions of . . . equilibrium, and macrostability, . . . which are inadequately understood. Especially on the present time scale of high-beta experiments, it is quite unlikely that microphenomena are important." Harold Grad (Director of the Magneto-Fluid Dynamics division of the Courant Institute of NYU), "Plasmas", Physics Today, December 1969, pp. 34-44.
(3) ". . . we define a `high density plasma` to mean a fully-ionized gas with density in excess of about 10.sup.16 particles/cm.sup.3 and with a temperature in excess of some 10.sup.16 .degree.K. Such a plasma is, generally speaking, collision-dominated and behaves in most respects as a fluid rather than as an assembly of independent particles. The instabilities to which the plasma is subject are consequently limited to the fluid, or magnetohydrodynamic instabilities; we are not concerned with the so-called velocity-space instabilities or those which derive from anisotropies in the pressure distribution." G. B. F. Niblett (of Culham Laboratory, UKAEA), "The Production and Containment of High Density Plasmas", Physics of Hot Plasmas, edited by Rye and Taylor, Plenum Press, 1970, pp. 361-378.
(4) "In a plasma the collision frequency is given by ##EQU104## and for
n.congruent.10.sup.16 /cm, T.congruent.10keV
this is of the order of 10.sup.6 /sec; hence for systems lasting about 1 sec, the distribution of particle velocities will be roughly isotropic and a magnetohydrodynamic description is appropriate." W. B. Thompson, "Introduction to Controlled Thermonuclear Research," pp. 259-271 of Plasma Physics (Reference [10.2]).
In the example just quoted, .nu.=0.64.times.10.sup.6 sec.sup.-1. To consider a hotter, denser plasma, postulate a confining field B of 20 tesla, density n=2.5.times.10.sup.16 cm.sup.-3 and temperature T=20keV in which case .nu.=0.56.times.10.sup.6 sec.sup.-1 and the magnetohydrodynamic description remains appropriate.
We shall employ the following state variables:
.rho.=fluid density
v=fluid velocity vector
.tau.=fluid temperature
B=magnetic field vector
E=electric field vector
j=electric current density of fluid
.rho..sub.e =electric charge density of fluid
n=normal vector to free fluid boundary
r=radius vector to a given position
We shall employ rationalized emu (electromagnetic units) together with the following parameters:
c=velocity of light in vacuum
c.sub.f =velocity of light in fluid=1/.sqroot..epsilon..mu..sub.e .congruent.c
c.sub.e =.mu..sub.e c.sub.f .congruent.c
.mu..sub.e =permeability of fluid (.mu..sub.e =1 in vacuum)
.epsilon.=dielectric constant of fluid
.rho..sub.e =(1/c.sub.3.sup.2).gradient..multidot.E
.eta..sub.1 =reciprocal of electron-ion collison frequency
.eta..sub.i =.eta..sub.i (.rho.,.tau.), (i=2,3,4,5) electrical coefficients proportional to .eta..sub.1
.eta.=.eta.(.rho.,.tau.)=electrical resistivity
.kappa.=.kappa.(.rho.,.tau.)=thermal conductivity
.mu.=.mu.(.rho.,.tau.)=hydrodynamic viscosity
.sigma..sub.H =source of heat, per unit volume
.sigma..sub.M =source of mass, per unit volume
.gamma.=c.sub.p /c.sub.V =ratio of specific heats (>2 for plasmas)
R=c.sub.p -c.sub.V =(.gamma.-1)c.sub.V >0=gas constant
p=R.rho..tau.=fluid pressure (perfect gas law)
We omit radiation pressure, although the present model may be readily modified to include it.
Propagation of Electromagnetism ##EQU105##
.gradient..multidot.B=0 (10.3)
Conservation of Charge ##EQU106##
Conservation of Mass ##EQU107##
Conservation of Momentum ##EQU108##
Conservation of Energy ##EQU109## where .gradient.v is the Jacobian matrix (.differential.v.sub.i /.differential.r.sub.j) of v with respect to the position vector r, and where * denotes matrix transposition.
Generalized Ohm's Law ##EQU110## where as usual ##EQU111##
Boundary Kinematics ##EQU112## which expresses the fact that the fluid-vacuum boundary surface moves with the fluid.
REFERENCES
[10.1] W. B. Thompson, Introduction to Plasma Physics, Addison-Wesley, 1962.
[10.2] IAEA, Plasma Physics, Vienna, 1965.
[10.3] L. Spitzer, Physics of Fully Ionized Gases, Interscience, 2nd ed., 1965.
[10.4] Shih-i Pai, Magneto-Gasdynamics and Plasma Dynamics, Prentice-Hall, 1962.
K. DYNAMIC STRUCTURAL STABILITY CRITERION
Introduction
The mathematical model of the plasma pseudo-equilibrium presented in sub-part J above will now be employed to analyze the dynamic stability of plasma configurations of the class disclosed in the present invention. It should be recalled that this analysis possesses assured validity only for dense, "collision-dominated" plasmas such as the medium beta and high beta plasmas involved in the preferred embodiments of this invention. For low beta plasmas the conducting fluid model is of questionable adequacy; for example, in a paper on low .beta. "toroidal diffuse pinch" configurations, Coppi and Rem [11.1], page 69) state that the continuum mechanical model "is no longer adequate at the high temperature typical of thermonuclear plasmas, where the mean free paths for interparticle collisons become longer than the length of the system." Fortunately, there are two independent considerations indicating that for useful power reactor purposes high beta is to be preferred; the choice of high beta then relegates the questionable adequacy of the present model for low beta to a topic of merely academic interest. The first of these reasons concerns dynamic stability, and the second concerns economic efficiency.
Despite two decades of heavily-funded intensive research by many nations, a feasible means of stable plasma confinement in a closed configuration (by means other than those of the present invention, first publicly mentioned [11.22]) has remained undiscovered. For example, Academician L. A. Arzimovich, director of Soviet fusion power research, asserts that ". . . stability of the plasma filament is . . . an unresolved problem" ([11.2], p. 1969). Likewise, H. R. Hulme, Chief of Nuclear Research for the British AEA, states that ". . . the most serious deficiency in our knowledge concerns the stability of the plasma and how it can be maintained" ([11.3], p. 110), and ". . . how to keep the plasma stable for a time long enough for thermonuclear reactions to take place is the outstanding problem" ([11.3], p. 92). Likewise another expert, D. A. Frank-Kamenetskii, states in 1972: "The main obstacle in the quest for the ideal magnetic trap is the instability problem" ([11.4], p. 4) and ". . . to achieve this goal one item is certainly necessary--this is finding a way to cope with plasma instabilities. This is the single difficulty, but it is so serious that no one can say how close the solution is or which way it lies." ([11.4], p. 6). [Emphases added in preceding quotations.]
It was shown below in sub-parts L and M that in order to design a plasma equilibrium which can be naturally stabilized by a high aspect ratio or artificially stabilized by feedback, one requires dynamic structural stability, which means dynamic behaviour that for small perturbations is relatively insensitive to variations in the system's specifications. (A plasma can have dynamic structural stability and yet be dynamically unstable, as will be explained below.) Accordingly, the most significant result of the present section will be a demonstration that, assuming static structural stability as a prerequisite, one always has dynamic structural stability to the extent that the continuum mechanical model used heretofor in sub-part J is validly applicable, i.e., to the extent that the plasma is of high beta and sufficiently dense to be collison-dominated.
It is well-known in this art (see e.g. [11.3], p. 97) that for a given fluid temperature and given confining magnetic field strength the thermonuclear power production per unit volume is proportional to .beta..sup.2. The volume of the magnetic field which must be generated is an important parameter of the cost of a fusion reactor. In terms of this particular cost-parameter, compare a low beta reactor of say .beta..sub.l =0.10 with a high beta reactor of say .beta..sub.h =0.90; since .beta..sub.h.sup.2 /.beta..sub.l.sup.2 =(0.9).sup.2 /(0.1).sup.2 =81, the low beta reactor would cost eighty-one times more to construct, for the same power output, as a high beta reactor. Of course, if for other reasons such as limitation of heat flux through the walls, additional cost of providing sufficient magnetic field and additional complexity to the mathematical model are acceptable in order to construct a low beta or medium beta reactor, the principles disclosed herein remain fully applicable to the construction of such a reactor with no essential changes.
Inadequacy of MFD Analyses
The present model will be called an EMTFD model. If one sets c=+.infin., i.e.
(1/c)=0 (11.1)
and ignores E, then one has an MTFD model. If one sets the transport coefficients equal to zero,
.eta.=.kappa.=.mu.=0 (11.2)
one obtains an MFD (Magneto-Fluid Dynamical) model known as the Lundquist equations. The MFD model has been the basis of virtually all prior stability analyses (see references [11.5] through [11.17]). One linearizes the MFD equations around an equilibrium solution, obtaining wave equations for normal modes of oscillation. Each normal mode evolves in time t proportionally to the real or imaginary part of
e.sup..omega.t, (0>t.fwdarw.+.infin.) (11.3)
where .omega. is the complex frequency of that mode. For MFD models .omega..sup.2 is necessarily a real number. If .omega..sup.2 >0, then there is one positive value of .omega. and the mode is unstable. If .omega..sup.2 >0, then .omega. is pure imaginary and this case is called "stable" in the prior literature. However, this terminology is unwarranted and the conclusions drawn therefrom are misleading (see FIG. 31). For example, U.S. Pat. No. 3,016,342 is based entirely upon this alleged "stability" criterion, as specified in its FIGS. 10 through 18, all of which are graphs of the criterion specified in its line 45 of page 31, but the validity of this theoretical criterion has been disproved by experimental tests. Indeed it is well-known in nonlinear mechanics [11.18] that unless .omega. has a non-zero real part, no valid conclusions regarding the true nonlinear system can be drawn from the linearized analysis. A simple example in this connection will prove instructive.
Simple Model of Anomalous Diffusion
According to accepted MFD "stability" theory, one would say that the equilibrium solution
x=0, y=0 (11.4)
of the elementary dynamical system.
x=y+(1/2).mu.(x.sup.2 +y.sup.2)x, (0>.mu.<1) (11.5a)
y=-x+(1/2).mu.(x.sup.2 +y.sup.2)y, (0>.mu.<1) (11.5b)
is "stable" because linearization around equilibrium yields
x=y, y=-x (11.6)
or, eliminating y,
x+x=0, x(t)=x.sub.o e.sup..omega.t, (11.7)
where ##EQU113## No matter how small .mu. may be, this "conclusion" is mathematically incorrect and quite misleading.
It is easy to prove [simply differentiate (11.9) and use the given dynamics (11.5a,b)] that ##EQU114## which means that x.sup.2 +y.sup.2 "blows up" to infinity faster than any exponential e.sup..gamma.t, .gamma.>0, no matter how large .gamma. may be chosen. In fact, the system blows up within a finite length of time, namely, at time: ##EQU115## This behavior is qualitatively different from anything that can be inferred from the usual linearized equations of physics or mechanics.
Note that, if the initial conditions are such that x.sub.o.sup.2 +y.sub.o.sup.2 is "small" relative to (1/.mu.), then the initial rate of motion away from equilibrium is very "slow." On the other hand, if such a "large" value as x.sub.o.sup.2 +y.sub.o.sup.2 =1/.mu. is a physically allowable initial ("perturbed") state, then the initial state will reach infinity in one second (truly this is a model of an "explosion").
Hence, in between the clearly allowable initial states and those too large to be considered, there must exist states with every possible growth rate between "small" and "large".
A study of (11.4) through (11.10) leads to the following conclusions:
(i) Linearized conservative mathematical models of plasma dynamics (i.e., MFD models) are intrinsically inadequate for prediction of observable stability.
(ii) Even though the supposedly negligible nonlinear terms are "infinitesimally small" for small initial perturbations, they permit an initially slowly-growing instability which could account for the phenomenon of "anomalous diffusion."
(iii) Neglect of the nonlinear terms means that it is impossible to compute the rate of the dispersion (the information necessary thereto, i.e., the value of .mu., has already been discarded).
Nonlinear Stability Theory [11.18]
Since James Clerk Maxwell's analysis of the stability of Watt's automatic governor a hundred years ago, the subject of stability of dynamical systems has been studied intensively by pure mathematicians and control system technologists. Important concepts have been introduced by such men as Routh (1877), Poincare (1885), Liapunov (1892), Hurwitz (1895), Malkin (1931), Persidskii (1933), Andronow and Pontrjagin (1937), Dubosin (1940), Krasovskii (1954), and Zubov (1957).
A logical terminology has been long established, and today in nonlinear mechanics and control theory it is universally accepted that practically useful results cannot be obtained without making the logical distinctions which are denoted by such terms as:
(Liapunov) stability (or "neutral stability");
local stability;
global stability;
asymptotic stability;
uniform asymptotic stability;
structural stability; and
total stability--also called
stability under continually acting disturbances.
In the present instance, one may classify the complex frequencies of (11.3) as follows:
______________________________________Re {.omega.} .ltoreq. 0 STABILITY (11.11a)Re {.omega.} > 0 INSTABILITY (11.11b)Re {.omega.} .noteq. 0 STRUCTURAL STABILITY (11.11c)Re {.omega.} < 0 ASYMPTOTIC STABILITY (11.11d)Re {.omega.} .ltoreq. - .delta. < 0 UNIFORM ASYMPTOTIC (11.11e) STABILITY______________________________________
where the word uniform implies that the degree of stability .delta. is independent of the mode under consideration.
Unless one is a priori assured that (11.11c) holds, i.e. unless the system is structurally stable, then the mere absence of instability, (11.11b), does not imply what is desired, namely (11.11e); in fact, the negation of (11.11c) would permit
Re{.omega.}=0 (11.12)
which, as shown in the preceding section, provides no guarantee whatsoever that true stability of the nonlinear system has been attained.
The difficulties associated with "stability" analyses based on (11.12) are depicted graphically in FIGS. 31, 32 and 33.
It is insufficient to attain (11.12) in the MFD model, as depicted in FIG. 31, because in the EMTFD model the complex frequencies on the imaginary axis will always move away from it as illustrated in FIG. 33. What is needed is a method of classifying the various regions of the imaginary axis as suggested in FIG. 32. This can be done by carrying out a complete EMTFD analysis and then letting the displacement current and the three transport coefficients tend to zero, as in (11.1) and (11.2), and noting whether the corresponding complex frequencies move onto the imaginary axis from the right half-plane or the left half-plane.
Linearization
It is an adequate approximation to assume
.kappa.=.kappa..sub.1 .tau..sup.5/2, .mu.=.mu..sub.1 .tau..sup.5/2, .eta.=.eta..sub.1 .tau..sup.-3/2 (11.13)
and to linearize around a uniform state
.rho.=.rho..sub.0, .tau.=.tau..sub.0, v=0, B=B.sup.0, E=0, j=0 (11.14)
where .rho..sub.0,.tau..sub.0 and B.sup.0 are constants. We employ cylindrical coordinates and take B.sup.0 parallel to the z-axis, i.e., B.sup.0 =.vertline.B.sup.0 .vertline.u.sub.z.
Define
.kappa..sub.0 =.kappa.(.tau..sub.0), .mu..sub.0 =.mu.(.tau..sub.0), .eta..sub.0 =.eta.(.tau..sub.0), (11.15)
and define kinematic viscosity .mu..sub.0 and thermal diffusivity .kappa..sub.0 by
.mu..sub.0 =.mu..sub.0 /.rho..sub.0, .kappa.=.kappa..sub.0 /.rho..sub.0 c.sub.p (11.16)
Then define perturbations (.rho.,.tau.,b,e) by
.rho.=.rho..sub.0 (1+.rho.), .tau.=.tau..sub.0 (1+.tau.), B=B.sup.0 +b, E=e, (11.17)
and define the Alfven speed .nu..sub.0 and the speed of sound .nu..sub.s by ##EQU116## Inserting the preceding quantities into the state equations in sub-part J one finds the following perturbation equations:
.differential..rho./.differential.t=-.gradient..multidot.v (11.19)
.differential..tau./.differential.t=-(.gamma.-1)(.gradient..multidot.v)+.gamma..chi..sub.0 .gradient..sup.2 .tau., (11.20) ##EQU117##
.differential.b/.differential.t=-.gradient..chi.e, (11.22)
(1/c.sup.2) (.differential.e/.differential.t)+j=.gradient..chi.b. (11.24)
where u.sub.z is the unit vector in the z-direction. (Here we have used a simplified Ohm's law based on the hypothesis that .eta..sub.1 =0, i.e. that the plasma is collision-dominated.)
Now consider a cylindrical plasma of arbitrary cross-section. Equations (11.19)-(11.24) can be solved in closed form for circular cross-sections by means of trigonometric and Bessel functions. If q denotes any of the state variables, assume the form
q=e.sup..omega.t+i(m.theta.+kz) q.sub.0 (r), (i.sup.2 =-1), (11.25)
and insert into (11.19)-(11.24). The result is the same as introduction of the operators ##EQU118## into the equations and manipulating them algebraically as if .omega.,.lambda..sup.2 and k were real or complex numbers. The most important waves correspond to eigenvalue problems of the form ##EQU119## Perturbation of the boundary conditions (1.5) and (1.10) determines perturbed boundary conditions which make (11.27) a well-posed problem defining Helmholtz numbers
.lambda..sup.2 =.lambda..sub.mk.sup.2 (11.28)
for every pair of integers (m,k), if one re-scales the z-axis so that the ends z=-.pi. and z=+.pi. are identified. For large aspect ratios this approximation of a torus by identifying the ends of a cylinder has adequate validity because in the present method the rotational transform is indeed zero and for large aspect ratios all field lines on the plasma surface are approximately straight lines parallel to the toroidal axis, and so any segment of the torus is adequately approximated by a segment of a plasma cylinder which is a portion of a straight theta pinch. The dispersion relation for a theta pinch in the MFD case is well-known (see equation (5-105) on page 135 of reference [11.11], as is the fact that according to MFD theory it is a "stable" configuration; moreover, a full EMTFD analysis yields precisely the same result, namely (1.11e), so that in this case the MFD prediction proves to be reliable.
The use of (11.26) in (11.19)-(11.24) leads to a polynomial dispersion relation of the form
.omega..sup.6 +.alpha..sub.5 .omega..sup.5 +.alpha..sub.4 .omega..sup.4 +.alpha..sub.3 .omega..sup.3 +.alpha..sub.2 .omega..sup.2 +.alpha..sub.1 .omega.+.alpha..sub.0 =0, (11.29)
where the coefficients .alpha..sub.i are polynomials in .mu..sup.2,k.sup.2,(1/c.sup.2), and the transport coefficients .eta..sub.0,.mu..sub.0,.kappa..sub.0. The six roots of (11.29) form three complex conjugate pairs corresponding to electromagnetosonic waves, "fast" magnetosonic waves, and "slow" magnetosonic waves (see [11.15] or [11.16]); equation (11.29) can be factored in closed form in the limit where (11.1) and (11.2) hold. In this case (11.12) holds if and only if
.lambda..sub.mk.sup.2 .gtoreq..phi.k.sup.2 (11.30a)
where by definition ##EQU120## In the second expression for .phi. we have made use of the boundary condition (1.10) in the readily-derived form ##EQU121## which enables one to plot .phi. as a function of .beta. as depicted in FIG. 7. In the MFD case the boundary conditions defining the Laplacian operator .gradient..sup.2 imply that it is self-adjoint, whence its eigenvalues .mu..sup.2 are necessarily real, and so (11.30) is automatically satisfied for all modes (m,k) whenever .phi.=0, i.e. whenever either of
.beta.=0, .beta.=1 (11.32)
holds. In the EMTFD case there is an analogue of (11.30a), in which (11.30b) is replaced by a more complicated function of .beta. and the transport coefficients, which behaves as illustrated in FIG. 8. The finite "jump" in .phi. at .beta.=0 is caused by the well-known "resistive instabilities", which constitute, as stated in paragraph two, page 1 of U.S. Pat. No. 3,278,384 the principal cause of the fact that the specific magnetic configurations provided by U.S. Pat. Nos. 2,910,414; 3,002,912, 3,012,955; 3,015,618; and 3,088,894 have all turned out to be dynamically unstable and therefore unworkable. However, in the EMTFD case .phi.=0 remains valid at
.beta.=1, (11.33)
as illustrated in FIG. 8.
Another indication of the desirability of high density comes from asymptotic analysis of the roots of (11.29) as k.sup.2 .fwdarw.+.infin.. Referring to (11.11e), this proves that ##EQU122## which shows (see Bass [11.23]) that improvements in the degree of asymptotic stability can, for short wave-length "kink" modes, be made only by increasing the ratio of the square of the density to the cube of the temperature.
As a check on the preceding analysis, note that if one neglects the displacement current, i.e. assumes (11.1), then the perturbation equations (11.19)-(11.24) become identical with those of Coppi [11.19] and Bernstein and Robinson [11.20].
The only author attempting to treat the case .eta..noteq.0, .kappa..noteq.0, .mu..noteq.0 seems to be Trehan [11.21], but his work is seriously flawed by several errors. In his equations (90)-(111), pages 37-39, he includes resistivity and viscosity; he then attempts to add thermal conductivity in his equations (123)-(131), pages 41-42, but all of these equations [except (124)] are incorrect because of omission of a factor .rho..sup..gamma. on the left-hand side of his equation (123). Correction of this error leads to major changes in his (127)-(128). (In reference [11.21], equations (127)-(128) are not solved correctly even as they stand, because comparison of (130)-(131) with the definition (99) implies that all plasma waves are growing exponentially.) When (129)-(130)-(131) of [11.21] are corrected, the results become special cases of the present results concerning (11.29).
In conclusion, the present analysis demonstrates that the class of magnetic confinement configurations disclosed in the present invention constitutes a class of high beta toroidal plasma equilibria which possess uniform asymptotic stability, against arbitrary small perturbations, provided only that the configurations are specified to have sufficiently large aspect ratios.
REFERENCES
[11.1] B. Coppi and J. Rem, Scientific American, vol. 227, July 1972, pp. 65-75.
[11.2] L. A. Arzimovich, Elementary Plasma Physics, Blaisdell, 1965 (translation from 1963 Russian edition).
[11.3] H. R. Hulme, Nuclear Fusion, Springer-Verlag New York, 1969.
[11.4] D. A. Frank-Kamenetskii, Plasma, the Fourth State of Matter, Plenum, 1972 (translation from 1968 Russian edition).
[11.5] W. B. Kunkel, ed. Plasma Physics in Theory and Application, McGraw-Hill, 1966.
[11.6] IAEA (Seminar Trieste, 1964), Plasma Physics, Vienna, 1965.
[11.7] L. Spitzer, Physics of Fully Ionized Gases, Interscience, 2nd ed. 1962.
[11.8] S. Glasstone and R. H. Lovberg, Controlled Thermonuclear Reactions, van Nostrand, 1960.
[11.9] D. J. Rose and M. Clark, Jr., Plasmas and Controlled Fusion, M.I.T. Press, 1961.
[11.10] W. B. Thompson, An Introduction to Plasma Physics, Addison-Wesley, 1962.
[11.11] G. Schmidt, Physics of High Temperature Plasmas, Academic Press, 1966.
[11.12] S. I. Pai, Magneto-Gasdynamics and Plasma Dynamics, Prentice-Hall, 1962.
[11.13] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Oxford University Press, 1961.
[11.14] W. I. Futterman, ed., Propagation and Instabilities in Plasmas, Stanford University Press, 1963.
[11.15] A. Jeffrey, Magnetohydrodynamics, Wiley, 1966.
[11.16] A. Jeffrey and T. Taniuti, Non-linear Wave Propagation, Academic Press, 1964.
[11.17] A. Jeffrey and T. Taniuti, Magnetohydrodynamic Stability and Thermonuclear Containment, Academic Press, 1966.
[11.18] W. Hahn, Stability of Motion, Springer-Verlag New York, 1967.
[11.19] B. Coppi, "On the Stability of Hydromagnetic Systems with Dissipation," reference [11.14], p. 70-86.
[11.20] I. B. Bernstein and B. B. Robinson, Ann. Phys., vol. 18 (1962), p. 110.
[11.21] S. K. Trehan, "Macroscopic Theory of Plasma Waves," pp. 21-46 of reference [11. ].
[11.22] R. W. Bass, paper 3A1, "The ULTRON Configuration," Bull. Amer. Phys. Soc., Nov. 1971, p. 1239.
[11.23] R. W. Bass, paper G3, "A Simple Necessary and Sufficient Condition for Structural Stability of Non-Conservative Toroidal Plasma Equilibria," Fifth Annual Meeting, Plasma Physics Div. of Amer. Phys. Soc., Nov. 1963.
L. MINIMUM ASPECT-RATIO FOR FINITE GYRO-RADIUS STABILIZATION
It was shown above in sub-part K that an infinitely long straight theta-pinch device possesses both static and dynamic structural stability. Hence it can be "bent" slightly, in a manner maintaining equilibrium, such as that disclosed in the present invention, without destruction of its dynamic stability. The objective of the present discussion is to calculate the maximum degree of such toroidal curvature which is tolerable.
It will be seen that for present purposes it suffices to approximate the plasma volume as an axially symmetric torus of circular cross-section. Let a denote the major radius and b denote the minor radius of this torus. The aspect ratio is defined to be the ratio (a/b).
Under physically reasonable assumptions, it will be demonstrated that the plasma equilibrium of FIGS. 2A, 2B, 3A and 3B possesses dynamic asymptotic stability if its aspect ratio exceeds approximately two hundred and forty.
For example, consider a confining field strength of B=20 tesla, at a plasma temperature of T=5keV, which corresponds to confinement of ionized hydrogen at a density of 10.sup.17 particles per cubic centimeter, and an average ion gyro-radius R.sub.c exceeding 0.05 cm, assuming a high beta plasma for a boundary layer or skin thickness L less than 0.10 cm, and a minor radius of b=1.0 cm, and adequate major radius would be a=2.4 m.
It is well-known in the present art that regions of unfavorable magnetic field curvature (concave toward the plasma) induce "local ripple instabilities" or "resistive ballooning instabilities." These "localized flute" modes are similar to that Rayleigh-Taylor instabilities present when a heavier liquid is supported by a lighter liquid; but in the latter case, such instabilities can be adequately damped by the presence of surface tension (Chandrasekhar [12.1]).
The electrostatic sheath surrounding the plasma has an effective surface tension (Rose and Clark [12.2], p. 291). It was shown by Rosenbluth, Krall and Rostoker [12.3] that this effect can damp local flute modes. Here the exposition and notation of Schmidt ([12.4], pp. 283-297) will be followed.
Without any significant loss of generality it will be assumed that a.gtoreq.(3/2L. Then equation (8-153) of Schmidt is violated and so his result regarding his equation (8-137) is applicable for 0<.beta..ltoreq.1. Note that in his notation R=a.
Assuming the field lines to be horizontal circles of radius of curvature R=a, the longest wavelength of a mode perpendicular to the field line would be about one-half of the minor circumference (because of the limit-cycles). Hence we may take, in Schmidt's notation, k=2.pi./.lambda. where .lambda..ltoreq..pi.b, or ##EQU123## It will further be assumed that ##EQU124##
Under these conditions, the hypothesis
a.gtoreq.240b (12.3)
implies that ##EQU125## which is precisely equation (8-137) of [12.4], and which therefore prevents the growth of local flute modes.
In order to prove (12.4), note that a.gtoreq.240b,R.sub.c .gtoreq.b/20 and (b/10).gtoreq.L together imply that ##EQU126## Now, using (12.5), and (12.1), one has ##EQU127## that is, (12.4) holds, as required.
REFERENCES
[12.1] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Oxford University Press, 1961.
[12.2] D. J. Rose and M. Clark, Jr., Plasmas and Controlled Fusion, M.I.T. Press, 1961.
[12.3] M. N. Rosenbluth, N. A. Krall and N. Rostoker, Nucl. Fusion, 1962 Suppl., Part 1, p. 143.
[12.4] G. Schmidt, Physics of High Temperature Plasmas, Academic Press, 1966.
M. SEVEN-STEP PROCEDURE FOR FEEDBACK STABILIZATION OF DISTRIBUTED-PARAMETER SYSTEMS SUCH AS PLASMAS
Sophisticated design of multi-channel stability augmentation or stabilization systems is not part of the prior art of the subject of the present invention. To date only very elementary applications of contemporary systems engineering theory to plasma stabilization have been reported. In such applications what constitutes negative feedback is obvious. However, when one considers simultaneously many degrees of freedom, as in a continuum-mechanical or distributed-parameter system, together with many sensors and many actuators, the problem of design of a feedback control filter network or equivalent digital computer-control means to link the sensors effectively to the actuators in such a multi-channel situation far exceeds the power of human intuition. It is then necessary to employ a systematic theoretical multi-channel design procedure, of the type known in the art of design and development of aerospace-vehicle stabilization systems, but hitherto little used in plasma confinement technology. The purpose of the present appendix is to indicate as concisely as possible how to apply such procedures to the problem of assuring adequate dynamic structural stability of plasmas confined by the principles of the present invention.
As shown in sub-parts E, F and G above, requiring static structural stability of the magnetic field line configuration on the plasma's surface leads to unfavorable curvature of the magnetic field lines on the outer part of that magnetic shell tangent to the plasma's surface. Thus the plasma is dynamically unstable (see sub-parts J and K), but the choice of static structural stability allows system identifiability (in a sense to be explained), which in turn allows the possibility of both observability and controllability with suitable choices of actuators and sensors; the latter properties then render it a routine matter to stabilize the naturally-unstable normal modes.
For each postulated kidney-bean shaped plasma cross-section 141, one may determine loci of conductors and corresponding strengths of currents in the conductors sufficient to produce a magneto-static equilibrium of the type defined here; analytical and numerical procedures for this determination have been presented in sub-parts C, D and I.
It was also shown in sub-parts J and K how to set up and linearize the EMTFD equations around an equilibrium and how to perform a normal mode analysis (small oscillation analysis) of it. This in effect converts the problem into the standard state-variable control problem (see Kalman, Falb, and Arbib [13.1]), except that there are infinitely many degrees of freedom. However, this problem may be "truncated" by neglect of all but a finite number of modes, say the first thirty modes. There are theoretical justifications for such a truncation (note the subject of completely continuous operators in the mathematical discipline of functional analysis) and it has been demonstrated experimentally to be effective both in computer simulations (see Dolph and Lewis ([13.2]) and in design of aerospace vehicles.
In terms of this high-dimensional but conventional control problem, one may now compute the matrices of controllability and observability defined by Kalman ([13.1]) and use them to select and size actuators and to select and position sensors and then to design a conventional feedback control system which will artificially stabilize the unstable modes (see for example ([13.3]).
In 1961 it was proved by R. W. Bass that if an open-loop autonomous linear system is both observable and controllable in the sense defined by R. E. Kalman in 1960, then constant coefficient gain factors can always be computed (see [13.1], p. 49) such that with their use physically realizeable passive stable networks or equivalent digital filters can be designed, by a computationally-effective algorithm, which will connect the sensors with the actuators in such a way that each and every pole (complex frequency of a natural mode) of the resulting closed-loop system can be made equal to any arbitrarily pre-specified value.
For high-beta plasmas, the local particle density may be measured by sensors based upon infrared lasers. Another standard sensor system consists of a plurality of magnetic probes (Rogowski loops or pick-up coils) situated in the vacuum surrounding the plasma and used to detect variations in the present vacuum field from the nominal or ideal field (that is, the field known theoretically to correspond with the equilibrium that it is desired to preserve). A corresponding system of actuators could be a plurality of secondary conductors in which secondary currents could be rapidly altered to generate alterations in the confining vacuum field; this is technologically feasible because there are available switches which can completely reverse the direction of a fifteen ampere current in as little as one hundred nanoseconds.
It has been disclosed in sub-parts J and K herein that the linearized EMTFD equations can be reduced in Helmholtz's equation (.gradient..sup.2 +.lambda..sup.2).phi.=0 for normal modes or standing waves of sound in an ideal gas. The complex normal frequencies are all roots of an associated polynomial dispersion relation of the tenth degree, whose coefficients are polynomials in the Helmholtz numbers (eigen-values of the Laplacian). The feedback loop is closed by noting that the Helmholtz numbers .lambda..sup.2 are determined by a boundary value problem involving the electric and magnetic field vectors on the external surface of the vacuum surrounding the plasma; but if these fields are switched in response to sensed variations in the field near the plasma, the eigenvalue problem for the Laplacian is no longer independent of the dispersion relation, but, together with it, constitutes an implicit system determining the complex frequencies. The feedback coefficients are free parameters in this system, which can be adjusted to give the complex frequencies desirable properties, for example negative real parts. One algorithm applicable for this purpose is the procedure of Kalman and Bass cited previously.
The technique outlined so far may be referred to as a prior model control. It is desirable to compare results obtained from such a priori models with those obtained from empirical models. Frequently the two approaches are pursued in parallel and ultimately merged, or else alternated until the results show sufficient similarity for an inference of reliability and adequacy to be obtained. Since the empirical model control is less well-known, it will be summarized here in greater detail in the remainder of part II in a form applicable to dynamic stabilization of statically structurally stable plasma equilibria.
Elements of Systems Engineering Theory:
a. Mathematical Modeling
Consider what engineers schematically and conceptually refer to as an arbitrary black box (fixed but unknown device or system).
Assume that information can be fed into and taken out of the box, but that otherwise the box is not influenced or communicated with.
A mathematical model is a rule for computation, or algorithm, which reliably relates the box's outputs to its inputs.
b. Empirical Model Diagnosis
Suppose that numerous input-output experiments are performed upon a given black box. Suppose that (with the aid of sophisticated mathematics and much numerical data-processing) one can deduce a satisfactory mathematical model of the input-output process represented by the box.
A model inferred from direct experiments (rather than deduced theoretically from detailed postulates about the specific character of the process inside the box) is called an empirical model. The procedure of performing input-output experiments and processing the data to identify the empirical model is called system diagnosis (see [13.8]). Contemporary plasma physics literature reflects little awareness of the revolutionary progress in empirical system diagnosis which has occurred in modern systems engineering.
c. Empirical Model Control
As already explained, a typical plasma equilibrium of interest lacks natural (dynamical) stability.
Suppose that one instruments the reactor with sensors which have the property that their dial pointers read zero if and only if the plasma is in the desired equilibrium. (For example, optical sensors can note displacements of the plasma volume; magnetic field variations can be detected by pick-up coils, etc.). The readings on these dials will constitute the system output.
New augment the system with actuators which can alter the forces acting on the plasma. (For example, additional strong currents can be switched on and off in coils near the plasma.)
The resultant experimental system now resembles an arbitrary black box. The system inputs are the deviations from zero of the sensor readings.
If an analog or digital computer control system, or computrol system, can be designed to link the sensors with the actuators via negative feedback, then the plasma can be forced to remain in or near the equilibrium, i.e. the equilibrium can be artificially stabilized. Specifically, if increasing outputs result in actuator commands which start the outputs to decrease again, then we have negative feedback. If such a computrol system can be designed and implemented, then the outputs will always fluctuate near zero, meaning that the plasma must remain near equilibrium; this is called automatic regulation.
The only difficulty in designing a computrol system is the following. With a single output to be regulated and only a single actuator to operate, it is easy to ascertain by direct physical intuition just exactly what constitutes negative feedback.
However, in a fusion reactor there may be many sensors and many actuators, and, what is worse, in the empirical model linking the uncontrolled output of the sensors to arbitrary actuator inputs, there may be scores of distinct, but cross-coupled channels. In such complicated systems, it is impossible to determine just what constitutes negative feedback by either intuition or by any classical engineering techniques.
Fortunately, higher mathematics and systems engineering theory, particularly under the stimulus of aerospace requirements, has in the past decade completely resolved this basic problem. Given the empirical system model, computer programs are available which will automatically unravel the above-mentioned complications and print out the complete design of the most effective computrol system to close the loop between sensor output and input by strictly negative feedback. When such a computrol system is developed, and used to link sensors with actuators, the information flow is circular and endless: a deviation of the system from equilibrium causes a sensor output which causes a computrol command to the actuator input which modifies the equilibrium, which thereby produces a new sensor output, etc. ad infinitum. Because of this endless circularity, one speaks of a closed loop system and says that one has developed a closed-loop automatic feedback control system, or feedback stabilization system.
In some cases the equilibrium of the process to be regulated has some natural stability, but not as much as would be desirable. In such cases, one speaks of a stability augmentation system rather than a stabilization system.
The preceding approach, called empirical model control, works well only to the extent that certain quantitative criteria are met. These criteria will now be explained.
d. Identifiability, Observability, Controllability
Obviously the preceding technical approach cannot be used if it is not possible to perform a satisfactory empirical diagnosis by means of input-output tests (accompanied by appropriate data processing), that is, if the given instrumented equilibrium lacks identifiability. It is easy to prove that an equilibrium lacks satisfactory identifiability if it lacks the property of static structural stability.
Once the parameters defining an empirical model have been identified by an experimental diagnostic procedure, it is possible to check whether or not two important quantitative properties, called observability and controllability are present in satisfactory amounts. Roughly speaking, the former measures the adequacy of the sensor subsystem and the latter does the same for the actuator subsystem. If numerical computation of either of these criteria indicates inadequacy, then the corresponding subsystem must be modified; subsequently, the open-loop system must be reidentified, and then the observability and controllability recomputed. This process must be repeated until adequate sensor and actuator subsystems have been developed. At that point the computrol algorithms already mentioned can be used to complete the design of an adequate computrol system.
When an engineering system has the property that small changes in the quantities defining the system can only produce small quantitative changes (and no qualitative changes) in the system's observable characteristics; then it has structural stability as defined above in sub-parts E, F and G.
A statically or dynamically structurally stable system can be in the previously-defined sense either stable or unstable; what structural stability means is that a mathematical model can be identified whose predictions are reliable even in the presence of unknown and undesired small disturbances, noise, manufacturing tolerances, operational errors, etc. Roughly speaking, a mathematical model has structural stability if it can be deformed slightly in arbitrary and unknown ways without altering the system's predicted behavior (be that either stability or instability).
For a magnetic confinement apparatus to be more than a delicate laboratory curiosity, it obviously must have structural stability (both static and dynamic).
Hitherto the most extensively analyzed magnetic confinement devices are the Stellarator and the Tokomak. Experiments have shown that neither has structural stability. The Stellarator at Garching bei Muenchen has exhibited a distribution of its radial density versus magnetic rotational transform angle which contradicts prior theoretical predictions; but it is easy to see that these predictions were based upon a model which fails to satisfy known criteria necessary for structural stability.
The same objection applies to the basic mathematical model of the Tokomak, which could not be made to work even briefly until, after a decade of trial-and-error efforts, it was learned that the confining magnetic field strength had to be within one hundredth of one per cent of its nominal design value in order for an equilibrium to be attained (see [13.4], p. 641, par. 3).
f. Experimental Evidence of Relevance of Static Structural Stability
The standard theory of toroidal magneto-fluid dynamical equilibria (which goes back to Spitzer, Kruskal, and Kulsrud [13.5] involves a nested family of toroidal magnetic shells, the most important characteristic of each shell being its rotational transform angle or Poincare rotation number [13.6]. When the number is rational at least one magnetic field line is closed; when it is irrational no field line is closed and every line winds ergodically over a dense subset of the entire shell. It can be proved theoretically that if the number is not constant and varies continuously with increases in the radial parameter indexing the shells, certain gross physical properties of the shells should jump, discontinuously, infinitely often as the radial parameter is increased. (A certain integral important in dynamic stability criteria may have value V.sub.1 at all rational values and value V.sub.2 at all irrational values.) Such paradoxical behavior contradicts all experience of macroscopic phenomena and indicates clearly the something is wrong with the mathematical model.
This suspicion of mathematical pathology has been strikingly confirmed by experimental results. On the basis of known results in those parts of electromechanical engineering studied by the discipline of nonlinear mechanics, one might expect a finite number of the rational shells to expand into shells of a finite thickness and to alternate with shells of finite thickness (consisting of shells of identical rotation numbers) having the properties associated with an irrational rotation number. (This is shown clearly in FIG. 8, page 43, [13.77].) To the surprise of many, this is exactly what has been observed in certain plasma physics experiments which were designed and performed with great care at Garching bei Muenchen (the Wendelstein l=2 Stellarator). As shown in FIG. 6, page 41, of [13.7], the plasma density drops sharply at alternating, finitely thick layers, centered on shells, the denominator of whose rotation number is 3, 4, 5, . . . , 15. This phenomena cannot be explained on the basis of the standard continuum-mechanical macroscopic theories. However, it would be very desirable to avoid the complexity of micro equilibrium studies and examinations of this topic based on kinetic theory.
g. Consequences for Fusion Reactor Configurations
Fortunately, the standard continuum-mechanical model of the preceding situation can be demonstrated to have glaring defects when it is examined from the point of view of the Andronov-Pontrjagin theory of structural stability. An important necessary (but far from sufficient) condition for the rotation number of a magnetic shell to remain invariant under small perturbations of the external currents generating the magnetic field, in which the values of the currents, or their geometry, or both, are altered sufficiently slightly, is that the rotation number must be rational. This shows two things. Firstly, there is to be expected a real physical distinction between the two types of shells. Secondly, none of the shells are components of a valid methematical model, with observable consequences, since the irrational shells all fail the rotation number test and the rational shells, though they pass this test, fail the further test that the only acceptable rotation number is identically zero.
The systems engineering approach thus automatically excludes every hitherto suggested closed toroidal magnetic confinement means.
h. Unique Acceptable Magnetic Confinement Means
There are two distinct aspects of the structural stability of a plasma-equilibrium model; these are, respectively, static and dynamic structural stability.
In the former, one poses a geometric or topological question, "Does the qualitative pattern of the magnetic field lines on the magnetic shell tangent to the plasma surface remain unchanged when small changes are made in the geometry of external current-carrying conductors, or in the values of the external currents?"
If not, it is unreasonable to expect to proceed reliably to the next question: "Given an equilibrium with fixed external conductors and currents, how does the plasma behave in time if its initial state is not exactly at, but merely near to, that of the equilibrium; and is this predicted time-behavior invariant under small changes in the external conductors and currents?"
As mentioned above in connection with the EMTFD analysis of sub-parts J and K, this second problem, called that of dynamic structural stability, can be attacked by techniques from the theory of small oscillations in fluid or continuous media. (Recall the standing acoustic modes in pipe organs or the standing microwave fields in cavity resonators.) Although plasma confinement researchers have exhibited no awareness of this decisively important fact, it is well known in non-linear stability theory that a linearized (or small oscillation) analysis can give accurate and reliable predictions of stability or instability regarding the true non-linear system only if the linearized model has dynamic structural stability. In other words, the prediction of whether or not a particular mode of oscillation will start to grow, or will instead decay in time, is valid for engineering purposes only if the linearized model has structural stability; but in sub-part K it was proved that the EMTFD model (unlike all prior models) does have this required property. [If one includes simultaneously displacement current, electrical resistivity, thermal conductivity, and viscosity in the one fluid EMTFD equations, and carries the Hall effect term and other finite gyro-radius terms in the complete version of Ohm's law due to Professor L. Spitzer, then it is a priori certain that the complex frequencies of all normal modes of oscillation will have non-zero real parts, which ensures dynamic structural stability.] Therefore, until recently the only remaining unsolved problem concerned static structural stability, the principal objective of the present invention.
A Seven-Step Procedure:
The following exposition is intended to give a brief review of how a modern system-theory and information-sciences oriented fusion reactor designer would go about controlling an unknown dynamical plant or process. The most important basic concepts are IDENTIFIABILITY, OBSERVABILITY, and CONTROLLABILITY.
Identifiability means that parameters defining an empirical input-output model of the process can be obtained from appropriate experimental test the testing generates data which can be processed by means of suitable computer programs of yield estimates of the parameters which specify the empirical mode. The requirement of static structural stability assures identifiability in the present method.
Observability and Controllability pertain to the selection, positioning, and sizing of sensors and actuators adequate for regulation of the plant.
Although particular sensing instruments and actuating means may be used in identifying the empirical model, they may be far from ideal with respect to effectiveness in designing a closed-loop regulator system for the process.
Observability and controllability are mathematical criteria due to R. E. Kalman [13.1] which enable the designer to assign quantitative figures-of-merit to various hypothetical choices of alternative sensor/actuator sub systems and thereby make an effective final choice.
Any sensor subsystem together with the empirical model of the plant constitutes a plant-sensor pair. A plant-actuator pair is defined similarly.
A plant-sensor pair is observable if at least one sensor is stimulated by every possible dynamical mode of the plant, i.e., if every variable defining the instantaneous state of the plant has an effect on at least one sensor. More specifically, the pair is observable if it is possible, in principle, to estimate with adequate effectiveness the complete set of state variables at any time by suitable computer-processing of the past sensor outputs (up to and including the present output). A plant sensor-pair is not observable if dynamical modes can be excited which have no effect upon any sensor, e.g., some internal state variable (such as a pressure) can be increasing exponentially yet no sensor will indicate that the plant is beginning to explode.
A plant-actuator pair is controllable if at least one actuator can excite every state variable. More specifically, the pair is controllable if it is possible, in principle, to force or drive any state variable to any desired value by appropriate command inputs to the actuators. A plant-actuator pair is not controllable if some dynamical modes exist upon which no actuator has any effect, e.g., the plant may be exploding but no possible combination of actuator inputs can halt or reverse the process.
After final selection of the sensor/actuator subsystems so that the plant-sensor pair has a high degree of observability and so that the plant-actuator pair has a high degree of controllability, it is possible to proceed with the design of a feedback control filter or control computer which will link the sensors to the actuators in a closed-loop feedback controlled system; specifically, the control computer embodies suitable information processing algorithms which enable the actuator inputs to respond to the sensor outputs in such a manner as to drive the plant's state toward its desired condition.
The unknown plant, process, or open-loop system may be identified by input-output experiments and appropriate data-processing in a sequence of three steps.
STEP ONE
In the first step, one identifies the Markov parameters of the unknown process. The details of two particular methods of doing this (a frequency-response method for processes which have natural input-output stability, and a regression-analysis type of method for naturally unstable processes) are given in [13.10].
STEP TWO
In the second step, one identifies the input-distribution-matrix, the output-distribution matrix, and the dynamical coefficient matrix, as well as the dimension n of the system's state-vector. These three matrices, together with n, can be computed from the Markov parameters by means of the Ho-Kalman Algorithm [13.1].
STEP THREE
In the third step, or check on step two, one obtains an independent estimate of the dimension n, and simultaneously greatly refined estimates of the precise numerical values of the three matrices, by means of techniques from multi-channel stochastic optimal nonlinear filtering [13.11]; compare [13.12] and [13.8].
At this point the diagnosis of the unknown plant may be regarded as complete.
STEP FOUR
The next step is to optimize the input- and the output-distribution matrices by means of imaginary experiments (performed analytically and computationally) in which the existing sensors and actuators are compared with other feasible sensor/actuator subsystems. For example, the nature of the sensors and actuators may be changed; or their location may be altered; or their associated gain parameters may be varied, and the results assessed quantitatively.
The quantitative measures of the sensor/actuator subsystem effectiveness are the magnitudes of the inverses of the matrices of controllability and observabiltiy. The first may be computed from the input-distribution matrix and the dynamical coefficient matrix, while the second may be computed from the output-distribution matrix and the dynamical coefficient matrix. A technique for optimizing the above-mentioned inverses, i.e., keeping the dynamical coefficient matrix fixed and optimizing the input- and output-distribution matrices, is given in [13.13].
It may turn out that the existing sensors and actuators are adequate. But if the preceding analysis indicates poor controllability or poor observability, or both, it will be necessary to modify the sensor/actuator subsystems and then repeat the first four steps above in order to make sure that the input-and output-distribution matrices have been properly identified and that the new system actually does have adequate controllability and observability.
STEP FIVE
The fifth step is to design (purely for the purpose of optimization studies, because it may not be practicable to fully implement it) an optimal feedback control filter based on optimal filtering theory and stochastic optimal control theory (cf. Mortensen [13.14] for the completely general nonlinear case). This step can be accomplished by means of the GASP computer program [13.15]; compare [13.3], Chapter 15, pp. 341-375.
STEP SIX
The sixth step is to design a practical version of the optimal feedback filter, in which the filter is a priori constrained as to the degree of complexity allowed.
The resultant feedback control filter may still not be practical; e.g., it may require impractically large gains (which result from poor controllability or poor observability). One cure is to improve the sensor/actuator subsystem, i.e., repeat steps five and six.
An alternative cure is to relax the constraint on the complexity of the control computer. Considering the revolutionary improvements in modern information-processing technology, this is often the technique of choice. It can be shown that if p is the number of poles in the control filter (which for economy it is desired to keep small), if l is the number of sensors, if m is the number of actuators and if n is the number of state variables (as determined in steps two and three), then step six can always be completed whenever ##EQU128##
In step five one finds the n optimal filter poles for a Kalman-Bucy filter. This is usually impractically large. In step six, if one is willing to take p as large as p=n-l then the two matrix equations become uncoupled and can always be solved, yielding a Luenberger filter. In the presently recommended approach one starts with the smallest integer p which is larger than [m/(m+l)](n-l) and performs step six. If simulations show that the resultant closed-loop system is not satisfactory, then replace p by (p+1) and repeat step six. When simulations reveal that the point of diminishing returns has been reached, i.e., further increases in the dimensionality of the filter yield comparatively insignificant performance improvements, then the value of p has been found which is best from the point of view of cost-effectiveness.
STEP SEVEN
The seventh and final step is to convert the system design into one which is still optimal when saturation is taken into account. The result is always a system in which some variable is switched in a discontinuous manner, i.e., a bang-coast-bang system in which a tri-stable switch can have either the positions +1, 0 or -1. Since bang-coast-bang systems have a name which connotes discontinuity, many non-specialists associate the name with an unacceptable kind of control, which makes it necessary to explain that in practice all closed-loop feedback control systems operate in a bang-coast-bang mode, when analyzed in sufficient detail. This is because some physical variable always saturates except when the closed-loop system's outputs are so close to their desired values that every component in the entire system is operating within its range of linearity, which rarely happens. The way around this difficulty is to take advantage of it.
Specifically, choose the m variables in the closed-loop system's actuator-subsystem's complete block-diagram which saturate first, and regard the saturation limits as a priori constraints on the size of the control inputs. In doing this it usually becomes necessary to raise the dimensionality n of the problem by incorporating the actual sensor subsystems and actual actuator subsystems together with the actual unknown plant or process in an experimental assembly known as an on-line subsystem test wherein the feedback control filter and those environmental interactions difficult to provide directly are provided by a hybrid computing system coupled directly with the other subsystems and operating in real time. Upon repeating steps one through six with this more detailed model, the dimension n is usually raised.
Now suppose that so far the input-output experiments upon the openloop plant have been conducted in an artificial laboratory environment not subject to various disturbances that would be encountered in practice. These effects should be incorporated (e.g., by computer simulation) in the on-linesubsystem test referred to above. Frequently, satisfactory results can be obtained by forcing only the unstable or least stable of the original n open-loop poles into new positions, and by forcing a few of the newly added p filter poles into desirable positions, and letting the remaining poles fall where they may.
It should be emphasized that by means of a separation loop it can be assured in advance that the filter poles are the same when the two subfilters comprising the feedback filter are considered separately as the open-loop components (to which steps one through three are applicable, e.g., for standardization purposes), and when the poles are regarded, after the loops are all closed, as part of the (n + p) poles of the entire closed-loop. (In conventional design procedures, they are always shifted, as there is no separation loop.)
The level of detail necessary to arrive at the bang-coast-bang model varies with the situation at hand. In the most detailed model, the saturating variable is the switching information itself (or output of the feedback control filter), and at this level there is almost perfect correspondence between the theoretical model and reality.
It was discovered by Madame Flugge-Lotz [13.16] that linear feedbacks to switches usually lead to the phenomenon of high-frequency switching called "chattering"; mathematically, the average of the values of each of the state variables behaves exactly as if the system were actually governed by perfectly proportional power amplifiers! The mathematical theory by which one can compute the equivalent linear system (corresponding to the actual switched chattering system) and vice versa was first derived in general by Andre and Seibert (see [13.17]and [13.18]). Mathematically, saturated control of this type forces the n state variables onto an (n-1) - dimensional subspace or "hyperplane" of the state space, and keeps the state on this hyperplane while "sliding" the state to its target value; for this reason Russian authors refer to "chattering " as operation in the sliding regime.
The foundations for a general technique for optimizing switching systems was laid by students of Lefschetz, especially Bellman, who pioneered in the variational treatment of general discontinuous control problems (cf. [13.19]). In 1956, R. W. Bass noted ([13.17], p. 191) that the optimal switcing times could be mechanized by an on-line simulation of the system "adjoint" to the original open-loop plant. In 1957, Pontriagin, Gamkrelidze, Boltianskii and Mischenko published the celebrated "Maximum Principle", which unifies and generalizes the previously known special cases into a complete theory [13.20]. The synthesis of optimal switching feedback systems was further developed by Wonham and Johnson [13.21] and by Bass, Gura and Webber [13.22], [13,23], [13.24].
The system as designed in step six must next be augmented by an additional fine-tuning filter, which processes the outputs of the sensors and the switching elements to give a correction feedback signal which modifies the output of the control feedback filter designed previously. (The use of feedback from switches as an input to the filters which drive the same switches is called "hard feedback" by Russian authors). The fine-tuning filter is a dual-mode filter (the modes being selected, by an additional switch, as high-frequency or low-frequency according to whether or not the main switches are chattering or stable, e.g., whether the valves are not saturated or are saturated).
One may now proceed to the completion of step seven by ultra-fine tuning. The actual switches will have the additional complexities of dead-zone, hysteresis, and time-delay. For techniques suitable for handling these effects, see Bass [13.17].
With the elimination or optimization of any remaining unacceptable steady-state "hunting", the seventh step is complete.
CONCLUSION: COMPUTER CONTROL
The seven steps outlined above will lead systematically to the design and development of a satisfactory computroller. Notice that the final version of the system which links the l sensor outputs with the m actuation inputs is a rather complicated collection of filters and switches. In principle each filter can be realized by a passive, stable RLC network, but active elements (summers and switches) are also needed. Thus one has both analog techniques and switches. The result begins to resemble a hybrid computer.
This suggests consideration of design of the computrol system, which links the sensor outputs (suitably transduced) to the actuator inputs, as a special-purpose digital filter or digital data-processing system. Another approach would be on-line use of a parallel-logic multi-processor general-purpose digital computer in real time. Considerations of cost-effectiveness relative to a particular application will determine this final decision.
REFERENCES
[13.1] R. E. Kalman, P. L. Falb, and M. Arbib, Topcis in Mathematical System Theory, McGraw Hill, 1969
[13.2] C. L. Dolph and D. C. Lewis, Quaterly of Applied Mathematics, vol. 16, 1958, pp. 97-110
[13.3] B. D. O. Anderson and J. B. Moore, Linear Optimal Control, Prentice-Hall, 1971
[13.4] New Scientist, Sept. 25, 1969
[13.5] L. Spitzer, Physics of Fully Ionized Gases, Interscience, 2nd ed. 1962
[13.6] E. A.Coddington and N. Levinson, Theory of Ordinary Differential Equations. McGraw Hill, 1955
[13.7] Physics Today, December 1969
[13.8] D. Graupe, Identification of Systems, van Nostrand-Reinhold, 1972
[13.9] W. Hahn, Stability of Motion, Springer-Verlag, New York 1967
[13.10] R. W. Bass, Diagnosis of Unknown Multichannel Processes from Input-Output Experiments, CSATR 72-031-01, Computer Software Analysts, Inc., Los Angeles, Calif. October 1972
[13.11] R. W. Bass, V. D. Norum and L. Schwartz, "Optimal Multichannel Nonlinear Filtering", J. Math. Analysis and Applications, vol. 16 , 1966, pp. 152-164
[13.12] R. S. Bucy and P. D. Joseph, Filtering for Stochastic Processes with Applications, Interscience, 1968
[13.13] C. D. Johnson, "Optimization of a Certain Quality of Complete Controllability and Observability for Linear Dynamical Systems", ASME Trans. J. Basic. Eng., vol. 91, 1961, pp. 228-238.
[13.14] R. E. Mortensen, "Stochastic Optimal Control With Noisy Observations", Int. J. of Control, vol. 4, 1966, pp. 455-464.
[13.15] R. F. Webber, W. Freested, and R. W. Bass, "The GASP Program--An Integrated Tool for Optimal Control and Filter Design", Proc. TACC, 1968.
[13.16] I. Flugge-Lotz, Discontinuous and Optimal Control, McGraw Hill, 1968
[13.17] R. W. Bass, "Equivalent linearization, nonlinear circuit synthesis and the stabilization and optimization of control systems, " Proceedings, Symposium on Nonlinear Circuit Analysis, Polytechnic Inst. of Brooklyn MRI Series vol. VI (1956) pp. 163-198
[13.18] J. Andre and P. Seibert, "Uber stuckweise lineare Differentialgleichungen, die bei Regelungsproblem en auftreten, I, II," Arch. Math., vol. 7, 1956, pp. 148-156, 157-165
[13.19] R. E. Bellman, I. V. Glickberg and O. A. Gross, Some Aspects of the Mathematical Theory of Control Processes, Rand Corp., Jan. 16, 1958
[13.20] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze and E. F. Mischenko, The Mathematical Theory of Optimal Processes, Interscience, 1962
[13.21] W. M. Wonham and C. D. Johnson, "Optimal Bang-Bang Control with Quadratic Performance Index," Trans. ASME Ser. D: J. Basic Eng., vol. 86 (1964), pp. 107-115
[13.22] R. W. Bass and I. Gura, "High Order System Design via State-Space Considerations," Proc. Joint Auto. Control Conf. (1965), pp. 311-318
[13.23] R. W. Bass and R. F. Webber, "On Synthesis of Optimal Bang-Bang Feedback Control Systems with Quadratic Performance Index," Proc. Joint Auto. Control Conf. (1965, pp. 213-219
[13.24] R. F. Webber and R. W. Bass, "Simplified Algebraic Characterization of Optimal Singular Control for Autonomons Linear Plants", Martin Marietta Corp. Res. Rep. R-67-6, March 1967
While the invention has been described in connection with specific embodiments thereof and in specific uses, various modifications thereof will occur to those skilled in the art without departing from the spirit and scope of the invention as set forth in the appended claims.
The terms and expressions which have been employed in this specification are used as terms of description and not of limitation, and there is no intention in the use of such terms and expressions to exclude any equivalents of the features shown and described or portions thereof, but it is recognized that various modifications are possible within the scope of the invention claimed.

Number	Name	Date
3009080	Loos	Nov 1961
3088894	Koenig	May 1963
3219534	Furth	Nov 1965
3278384	Lenard et al.	Oct 1966
3663362	Stix	May 1972
3677890	Hartman	Jul 1972
3692626	Ohkawa	Sep 1972

Confinement of high temperature plasmas

Information

Patent Number

Date Filed

Date Issued

Inventors

Original Assignees

Examiners

Agents

CPC

US Classifications

Field of Search

US

International Classifications

Abstract

Description

Claims

BACKGROUND OF THE INVENTION

US Referenced Citations (7)

Non-Patent Literature Citations (7)

Continuation in Parts (1)

Entry
Soviet Physics-Technical Physics, vol. 7, No. 3, 9/62, pp. 187-191, by Skornyakov.
Physics of Fluids, vol. 9, 1966, pp. 2295-2296, by Yoshikawa et al.
Technology Review 12/76, pp. 20-24, 33-39, 41-43.
MATT-1050, 8/74, pp. 526-529.
ERDA-28, 1/75, pp. 1-3, 8-10.
Bulletin of the American Physical Society, 11/16/71, p. 1239, by Bass.
Proceedings of the Utah Academy of Science, Arts & Letters, vol. 50, part 2, 1973, by Gardner et al., pp. 1-11.