ISMFA - International School on Magnetohydrodynamics and Fusion Applications

09-16 September 2011

---

Plasma Focus Numerical Experiments

knowledge should be freely accessible to all

Sor Heoh Saw and Sing Lee in collaboration with Erol Kurt

Reference to the Lee model code should be given as follows:

Lee S. Radiative Dense Plasma Focus Computation Package (2011): RADPF www.plasmafocus.net http://www.intimal.edu.my/school/fas/UFLF/

Module 1: Introduction - The Plasma Focus and the Lee model

Summary

1.1 Description of the plasma focus; how it works, dimensions and lifetimes of the hot dense plasma

1.2 Scaling properties of the plasma focus

1.3 The radiative Lee model: the 5 phases

1.4 Using the Lee model as reference for diagnostics

1.5 Insights on plasma focus from numerical experiments using Lee model code

Description of the plasma focus. How it works, dimensions and lifetimes of the hot dense plasma

1.1.1 Introduction

The Plasma Focus is a compact powerful pulsed source of multi-radiation [1]. Even a small table- top sized 3 kJ plasma focus produces an intense burst of radiation with extremely high powers. For example when operated in neon, the x-ray emission power peaks at 10⁹ W over a period of nanoseconds. When operated in deuterium the fusion neutron burst produces rates of neutron typically 10¹⁵ neutrons per second over burst durations of tens of nanosecond. The emission comes from a point source making these devices among the most powerful laboratory pulsed radiation sources in the world. These sources are plasma-based.

When matter is heated to a high enough temperature, it ionizes and becomes plasma. It emits electromagnetic radiation. The spectrum depends on the temperature and the material. The higher the temperature and the denser the matter, the more intense is the radiation. Beams of electrons and ions may also be emitted. If the material is deuterium, nuclear fusion may take place if the density and temperature are high enough. In that case neutrons are also emitted. Typically the temperatures are above several million K and compressed densities above atmospheric density starting with a gas a hundredth of an atmospheric density.

One way of achieving such highly heated material is by means of an electrical discharge through gases. As the gas is heated, it expands, lowering the density and making it difficult to heat further. Thus it is necessary to compress the gas whilst heating it, in order to achieve sufficiently intense conditions. An electrical discharge between two electrodes produces an azimuthal magnetic field which interacts with the column of current, giving rise to a self compression force which tends to constrict (or pinch) the column. In order to ‘pinch’, or hold together, a column of gas at about atmospheric density at a temperature of 1 million K, a rather large pressure has to be exerted by the pinching magnetic field. Thus an electric current of at least hundreds of kA are required even for a column of small radius of say 1 mm. Moreover the dynamic process requires that the current rises very rapidly, typically in under 0.1 ms in order to have a sufficiently hot and dense pinch. Such a pinch is known as a super-fast super-dense pinch; and requires special MA fast-rise (ns) pulsed-lines. These lines may be powered by capacitor banks, and suffer the disadvantage of conversion losses and high cost due to the cost of the high technology pulse-shaping line, in addition to the capacitor banks.

A superior method of producing the super-dense and super-hot pinch is to use the plasma focus. Not only does this device produce superior densities and temperatures, moreover its method of operation does away with the extra layer of technology required by the expensive and inefficient pulse-shaping line. A simple capacitor discharge is sufficient to power the plasma focus.

1.1.2 The plasma focus

The plasma focus is divided into two sections. The first is a pre-pinch (axial) section. The function of this section is primarily to delay the pinch until the capacitor discharge (rising in a damped sinusoidal fashion) approaches its maximum current. This is done by driving a current sheet down an axial (acceleration) section until the capacitor current approaches its peak. Then the current sheet is allowed to undergo transition into a radial compression phase. Thus the pinch starts and occurs at the top of the current pulse. This is equivalent to driving the pinch with a super-fast rising current; without necessitating the fast line technology. Moreover the intensity which is achieved is superior to the line driven pinch.

Figure 1. Schematic of the axial and radial phases. The left section depicts the axial phase, the right section the radial phase. In the left section, z is the effective position of the current sheath-shock front structure. In the right section r_s is the position of the inward moving shock front driven by the piston at position r_p. Between r_s and r_p is the radially imploding slug, elongating with a length z_f. The capacitor, static inductance and switch powering the plasma focus is shown for the axial phase schematic only.

The two-phase mechanism of the plasma focus [1] is shown in figure 1. The inner electrode (anode) is separated from the outer concentric cathode by an insulating backwall. The electrodes are enclosed in a chamber, evacuated and typically filled with gas at about 1/100 of atmospheric pressure. When the capacitor voltage is switched onto the focus tube, breakdown occurs axisymmetrically between the anode and cathode across the backwall. The ‘sheet’ of current lifts off the backwall as the magnetic field (B_q) and it’s inducing current (J_r) rises to a sufficient value.

Axial phase: The J_r x B_q force then pushes the current sheet, accelerating it supersonically down the tube. This is very similar to the mechanism of a linear motor. The speed of the current sheet, the length of the tube and the rise time of the capacitor discharge are matched so that the current sheet reaches the end of the axial section just as the discharge reaches its quarter cycle. This phase typically lasts 1-3 ms for a plasma focus of several kJ.

Radial Phase: The part of the current sheet in sliding contact with the anode then ‘slips’ off the end ‘face’ of the anode forming a cylinder of current, which is then pinched inwards. The wall of the imploding plasma cylinder has two boundaries (see figure 1 radial phase). The inner face of the wall, of radius r_s is an imploding shock front. The outer side of the wall, of radius r_p is the imploding current sheet, or magnetic piston. Between the shock front and the magnetic piston is the annular layer of plasma. Imploding inwards at higher and higher speeds, the shock front coalesces on-axis and a super-dense, super-hot plasma column is pinched onto the axis (see figure 2 [2]). This column stays super-hot and super-dense for typically ten ns for a small focus. The column then breaks up and explodes. For a small plasma focus of several kJ, the most intense emission phase lasts for the order of several ns. The radiation source is spot-like (1mm diameter) when viewed end-on.

Figure 2. Dense plasma focus device. Image from Glenn Millam. Source: Focus Fusion Society For an animation of the plasma focus click here.

1.1.3 Radial dynamics of the plasma focus

Figure 3 shows a drawing of a typical plasma focus, powered by a single capacitor, switched by a simple parallel-plate spark-gap. The anode may be a hollow copper tube so that during the radial pinching phase the plasma not only elongates away from the anode face but also extends and elongates into the hollow anode (see figure 2). In figure 3 is shown the section where the current sheet is accelerated axially and also the radial section. Also shown in the same figure are shadowgraphs [3] taken of the actual radially imploding current sheet-shock front structure. The shadowgraphs are taken in a sequence, at different times. The times indicated on the shadowgraphs are relative to the moment judged to be the moment of maximum compression.

That moment is taken as t=0. The quality of the plasma compression can be seen to be very good, with excellent axisymmetry, and a very well compressed dense phase. In the lower left of figure 3 are shown the current and voltage signatures of the radial implosion [4], occurring at peak current. The implosion speeds are measured and has a peak value approaching 30 cm/ms.

This agrees with modelling, and by considering shock wave theory together with modelling [5] of subsequent reflected shock wave and compressive effects, a temperature of 6 million K (0.5 keV) is estimated for the column at peak compression, with a density of 2x10¹⁹ ions per cm³. The values quoted here are for the UNU/ICTP PFF 3 kJ device.

Figure 3. UNU/ICTP PFF - Design, Signatures and Dynamics

1.2 Scaling properties of the plasma focus

1.2.1 Various plasma focus devices

In figure 4a is shown the UNU ICTP PFF 3 kJ device [4-6] mounted on a 1 m by 1 m by 0.5 m trolley, which was wheeled around the International Centre for Theoretical Physics (ICTP) for the 1991 and 1993 Plasma Physics Colleges during the experimental sessions. The single capacitor is seen in the picture mounted on the trolley. In contrast, figure 4b shows the PF1000, the 1 MJ device [7] at the International Centre for Dense Magnetised Plasmas (ICDMP) in Warsaw, Poland. Only the chamber and the cables connecting the plasma focus to the capacitors are shown. The capacitor bank with its 288 capacitors, switches and chargers are located in a separate hall.

Figure 4a. 3kJ UNU ICTP PFF Figure 4b. 1 MJ PF1000 plasma focus

We show here the characteristics of several plasma focus devices [7].

Table 1 Characteristics of three plasma focus devices

Plasma Focus Devices	E0 (kJ)	a (cm)	Z0 (cm)	V0 (kV)	P0 (Torr)	Ipeak (kA)	va (cm/us)	ID (kA/cm)	SF [(kAcm^-1) Torr^0.5]	Yn (10 ^8)
PF1000	486	11.6	60	27	4	1850	11	160	85	1100
UNU ICTP PFF	2.7	1.0	15.5	14	3	164	9	173	100	0.20
PF400J	0.4	0.6	1.7	28	7	126	9	210	82	0.01

In table 1 we look at the PF1000 and study its properties at typical operation with device storage at 500 kJ level. We compare this big focus with two small devices at the kJ level.

From table 1 we note:

Voltage and pressure do not have any particular relationship to E₀.

Peak current I_peak increases with E₀.

Anode radius ‘a’ increases with E₀.

ID (current per cm of anode radius) I_peak/a is in a narrow range from 160 to 210 kA/cm

SF (speed or drive factor) (I_peak/a)/P₀^0.5 is 82 to 100 kAcm^-1/Torr^0.5 deuterium gas [8].

Peak axial speed v_a is in the narrow range 9 to 11 cm/us.

Fusion neutron yield Y_n ranges from 10⁶ for the smallest device to 10¹¹ for the PF1000.

We stress that whereas the ID and SF are practically constant at around 180 kA/cm and (90 kA/cm)/Torr^0.5 deuterium gas throughout the range of small to big devices, Y_n changes over 5 orders of magnitude.

The data of table 1 is generated from numerical experiments [5,9] and most of the data has been confirmed by actual experimental measurements and observations.

Table 2 Properties of three plasma focus devices

Plasma focus Devices	c= b/a	a (cm)	Tpinch (10^6/K)	vp (cm/ms)	rmin (cm)	zmax (cm)	Pinch Duration (ns)	rmin/a	zmax/a	Pinch duration/a (ns/cm)
PF1000	1.4	11.6	2	13	2.2	19	165	0.17	1.6	14
UNU ICTP PFF	3.4	1.0	8	26	0.13	1.4	7.3	0.14	1.4	8
PF400J	2.6	0.6	6	23	0.09	0.8	5.2	0.14	1.4	9

Table 2 compares the properties of a range of plasma focus devices. The properties being compared in this table are mainly related to the radial phase.

From table 2 we note:

i. The pinch temperature T_pinch is strongly correlated to the square of the radial pinch speed v_p.

ii. The radial pinch speed v_p itself is closely correlated to the value of v_a and c=b/a; so that for a constant v_a, v_p is almost proportional to the value of c.

iii. The dimensions and lifetime of the focus pinch scale as the anode radius ‘a’.

r_min/a (almost constant at 0.14-0.17)

z_max/a (almost constant at 1.5)

iv. Pinch duration has a relatively narrow range of 8-14 ns per cm of anode radius.

v. The pinch duration per unit anode radius is correlated to the inverse of T_pinch.

T_pinch itself is a measure of the energy per unit mass. It is quite remarkable that this energy density at the focus pinch varies so little (factor of 5) over a range of device energy of more than 3 orders of magnitude.

This practically constant pinch energy density (per unit mass) is related to the constancy of the axial speed moderated by the effect of the values of c on the radial speed.

The constancy of r_min/a suggests that the devices also produce the same compression of ambient density to maximum pinch density; with the ratio (maximum pinch density)/ (ambient density) being proportional to (a/r_min)². So for two devices of different sizes starting with the same ambient fill density, the maximum pinch density would be the same.

From the above discussion, we may put down as rule-of-thumb the following scaling relationships, subject to minor variations caused primarily by the variation in c.

i. Axial phase energy density (per unit mass) constant

ii. Radial phase energy density (per unit mass) constant

iii. Pinch radius ratio constant

iv. Pinch length ratio constant

v. Pinch duration per unit anode radius constant

1.2.2 Summarising

i. The dense hot plasma pinch of a small E₀ plasma focus and that of a big E₀ plasma focus have essentially the same energy density, and the same mass density.

ii. The big E₀ plasma focus has a bigger physical size and a bigger discharge current. The size of the plasma pinch scales proportionately to the current and to the anode radius, as does the duration of the plasma pinch.

iii. The bigger E₀, the bigger ‘a’, the bigger I_peak, the larger the plasma pinch and the longer the duration of the plasma pinch. The larger size and longer duration of the big E₀ plasma pinch are essentially the properties leading to the bigger neutron yield compared to the yield of the small E₀ plasma focus.

The above description of the plasma focus combines data from numerical experiments, consistent with laboratory observations.

The next section describes the Lee model code.

1.3 The radiative Lee model: the 5 phases

The Lee model couples the electrical circuit with plasma focus dynamics, thermodynamics, and radiation, enabling a realistic simulation of all gross focus properties. The basic model, described in 1984 [1], was successfully used to assist several projects [4-6]. Radiation-coupled dynamics was included in the five-phase code, leading to numerical experiments on radiation cooling [5]. The vital role of a finite small disturbance speed discussed by Potter in a Z-pinch situation [10] was incorporated together with real gas thermodynamics and radiation-yield terms. This version of the code assisted other research projects [4,8,11,12] and was web published in 2000 [13] and 2005 [14]. Plasma self-absorption was included in 2007 [13], improving the SXR yield simulation. The code has been used extensively in several machines including UNU/ICTP PFF [3,8,11,12], NX2 [12,15-17], and NX1 [15,18] and has been adapted for the Filippov-type plasma focus DENA [19]. A recent development is the inclusion of the neutron yield Y_n using a beam–target mechanism [20-24], incorporated in recent versions [5] of the code (versions later than RADPFV5.13), resulting in realistic Y_n scaling with I_pinch [20,21]. The versatility and utility of the model are demonstrated in its clear distinction of I_pinch from I_peak [25] and the recent uncovering of a plasma focus pinch current limitation effect [22,23], as static inductance is reduced towards zero. Extensive numerical experiments had been carried out systematically resulting in the uncovering of neutron [20,21] and SXR [26-33] scaling laws over a wider range of energies and currents than attempted before. The numerical experiments also gave insight into the nature and cause of ‘neutron saturation [9,30,34]. The description, theory, code, and a broad range of results of this “Universal Plasma Focus Laboratory Facility” are available for download from [5].

A brief description of the 5-phase Lee model is given in the following.

1.3.1 The 5 phases

The five phases (a-e) are summarised [5,13,14,27, 31,33,35] as follows:

a. Axial Phase (see figure 1 left part)

Described by a snowplow model with an equation of motion which is coupled to a circuit equation. The equation of motion incorporates the axial phase model parameters: mass and current factors f_m and f_c. The mass swept-up factor f_m accounts for not only the porosity of the current sheet but also for the inclination of the moving current sheet-shock front structure, boundary layer effects, and all other unspecified effects which have effects equivalent to increasing or reducing the amount of mass in the moving structure, during the axial phase. The current factor f_c accounts for the fraction of current effectively flowing in the moving structure (due to all effects such as current shedding at or near the back-wall, and current sheet inclination). This defines the fraction of current effectively driving the structure, during the axial phase.

Figure 5. Schematic of radius vs time trajectories to illustrate the radial inward shock phase when r_s moves radially inwards, the reflected shock (RS) phase when the reflected shock moves radially outwards, until it hits the incoming piston r_p leading to the start of the pinch phase (t_f) and finally the expanded column phase.

b. Radial Inward Shock Phase (see figure 1 right part, also figure 2)

Described by four coupled equations using an elongating slug model. The first equation computes the radial inward shock speed from the driving magnetic pressure. The second equation computes the axial elongation speed of the column. The third equation computes the speed of the current sheath, (magnetic piston), allowing the current sheath to separate from the shock front by applying an adiabatic approximation [5,7]. The fourth is the circuit equation. Thermodynamic effects due to ionization and excitation are incorporated into these equations, these effects being particularly important for gases other than hydrogen and deuterium. Temperature and number densities are computed during this phase using shock-jump equations. A communication delay between shock front and current sheath due to the finite small disturbance speed [10,35] is crucially implemented in this phase. The model parameters, radial phase mass swept-up and current factors f_mr and f_cr are incorporated in all three radial phases. The mass swept-up factor f_mr accounts for all mechanisms which have effects equivalent to increasing or reducing the amount of mass in the moving slug, during the radial phase. The current factor f_cr accounts for the fraction of current effectively flowing in the moving piston forming the back of the slug (due to all effects). This defines the fraction of current effectively driving the radial slug.

c. Radial Reflected Shock (RS) Phase (See figure 5)

When the shock front hits the axis, because the focus plasma is collisional, a reflected shock develops which moves radially outwards, whilst the radial current sheath piston continues to move inwards. Four coupled equations are also used to describe this phase, these being for the reflected shock moving radially outwards, the piston moving radially inwards, the elongation of the annular column and the circuit. The same model parameters f_mr and f_cr are used as in the previous radial phase. The plasma temperature behind the reflected shock undergoes a jump by a factor close to 2. Number densities are also computed using the reflected shock jump equations.

d. Slow Compression (Quiescent) or Pinch Phase (See figure 5)

When the out-going reflected shock hits the inward moving piston, the compression enters a radiative phase in which for gases such as neon, radiation emission may actually enhance the compression where we have included energy loss/gain terms from Joule heating and radiation losses into the piston equation of motion. Three coupled equations describe this phase; these being the piston radial motion equation, the pinch column elongation equation and the circuit equation, incorporating the same model parameters as in the previous two phases. The duration of this slow compression phase is set as the time of transit of small disturbances across the pinched plasma column. The computation of this phase is terminated at the end of this duration.

e. Expanded Column Phase

To simulate the current trace beyond this point we allow the column to suddenly attain the radius of the anode, and use the expanded column inductance for further integration. In this final phase the snow plow model is used, and two coupled equations are used similar to the axial phase above. This phase is not considered important as it occurs after the focus pinch.

We note [31] that in radial phases b, c and d, axial acceleration and ejection of mass caused by necking curvatures of the pinching current sheath result in time-dependent strongly center-peaked density distributions. Moreover the transition from phase d to phase e is observed in laboratory measurements to occur in an extremely short time with plasma/current disruptions resulting in localized regions of high densities and temperatures. These centre-peaking density effects and localized regions are not modeled in the code, which consequently computes only an average uniform density, and an average uniform temperature which are considerably lower than measured peak density and temperature. However, because the four model parameters are obtained by fitting the computed total current waveform to the measured total current waveform, the model incorporates the energy and mass balances equivalent, at least in the gross sense, to all the processes which are not even specifically modeled. Hence the computed gross features such as speeds and trajectories and integrated soft x-ray yields have been extensively tested in numerical experiments for several machines and are found to be comparable with measured values.

1.4 Using the Lee model as reference for diagnostics

1.4.1 From measured current waveform to modelling for diagnostics

The Lee model code [5,13,14] is configured [9] to work as any plasma focus by inputting:

Bank parameters, L₀, C₀ and stray circuit resistance r₀;

Tube parameters b, a and z₀;

Operational parameters V₀ and P₀ and the fill gas.

The computed total current waveform is fitted to the measured waveform by varying model parameters f_m, f_c, f_mr and f_cr one by one, until the computed waveform agrees with the measured waveform.

First, the axial model factors f_m, f_c are adjusted (fitted) until the features in figure 6: ‘1’ computed rising slope of the total current trace; ‘2’ the rounding off of the peak current as well as ‘3’ the peak current itself are in reasonable (typically very good) fit with the measured total current trace (see figure 6, measured trace fitted with computed trace).

Then we proceed to adjust (fit) the radial phase model factors f_mr and f_cr until features ‘4’ the computed slope and ‘5’ the depth of the dip agree with the measured values. Note that the fitting of the computed trace with the measured current trace is done up to the end of the radial phase which is typically at the bottom of the current dip. Fitting of the computed and measured current traces beyond this point is not done. If there is significant divergence of the computed with the measured trace beyond the end of the radial phase, this divergence is not considered important.

In this case, after fitting the five features ‘1’ to ‘5’ above, the following fitted model parameters are obtained: f_m=0.1, f_c=0.7, f_mr=0.12, f_cr=0.68.

From experience it is known that the current trace of the focus is one of the best indicators of gross performance. The axial and radial phase dynamics and the crucial energy transfer into the focus pinch are among the important information that is quickly apparent from the current trace [27-29].

Figure 6. The 5-point fitting of computed current trace to measured (reference) current trace. Point 1 is the current rise slope. Point 2 is the topping profile. Point 3 is the peak value of the current. Point 4 is the slope of the current dip. Point 5 is the bottom of the current dip. Fitting is done up to point 5 only. Further agreement or divergence of the computed trace with/from the measured trace is only incidental and not considered to be important.

The exact time profile of the total current trace is governed by the bank parameters, by the focus tube geometry and the operational parameters. It also depends on the fraction of mass swept-up and the fraction of sheath current and the variation of these fractions through the axial and radial phases. These parameters determine the axial and radial dynamics, specifically the axial and radial speeds which in turn affect the profile and magnitudes of the discharge current.

There are many underlying mechanisms in the axial phase such as shock front and current sheet structure, porosity and inclination, boundary layer effects and current shunting and fragmenting which are not simply modeled; likewise in the radial phase mechanisms such as current sheet curvatures and necking leading to axial acceleration and ejection of mass, and plasma/current disruptions. These effects may give rise to localized regions of high density and temperatures. The detailed profile of the discharge current is influenced by these effects and during the pinch phase also reflects the Joule heating and radiative yields. At the end of the pinch phase the total current profile also reflects the sudden transition of the current flow from a constricted pinch to a large column flow. Thus the discharge current powers all dynamic, electrodynamic, thermodynamic and radiation processes in the various phases of the plasma focus. Conversely all the dynamic, electrodynamic, thermodynamic and radiation processes in the various phases of the plasma focus affect the discharge current. It is then no exaggeration to say that the discharge current waveform contains information on all the dynamic, electrodynamic, thermodynamic and radiation processes that occur in the various phases of the plasma focus. This explains the importance attached to matching the computed total current trace to the measured total current trace in the procedure adopted by the Lee model code. Once matched, the fitted model parameters assure that the computation proceeds with all physical mechanisms accounted for, at least in the gross energy and mass balance sense.

1.4.2 Diagnostics-time histories of dynamics, energies and plasma properties computed from the measured total current waveform by the code

During every adjustment of each of the model parameters the code goes through the whole cycle of computation. In the last adjustment, when the computed total current trace is judged to be reasonably well fitted in all 5 waveform features, computed time histories are presented, in figure 7a-7o as an example, as follows: for the NX2 operated at 11 kV, 2.6 Torr neon [9,33].

Figure 7a. Fitted computed I_total Figure 7b. Computed I_total & I_plasma

Figure 7c. Tube voltage Figure 7d. Axial trajectory and speed

Figure 7e. Radial trajectories Figure 7f. Length of elongating structure

Figure 7g. Speeds in radial phases Figure 7h. Tube inductance-axial & radial phases

Figure 7i. Total inductive energy Figure 7j. Piston work and DR energy; both traces overlap

Figure 7k. DR axial and radial phases Figure 7l. Peak & averaged uniform n_i

Figure 7m. Peak & averaged uniform n_e Figure 7n.Peak and averaged uniform T

Figure 7o. Neon Soft x-ray power

1.4.3 Comments on computed quantities by Lee model code

i. The computed total current trace typically agrees very well with the measured because of the fitting. The end of the radial phase is indicated in Figure 7a. Plasma currents are rarely measured. We had done a comparison of the computed plasma current with measured plasma current for the Stuttgart PF78 which shows good agreement of our computed to the measured plasma current [28]. The computed plasma current in this case of the NX2 is shown in figure 7b.

ii. The computed tube voltage is difficult to compare with measured tube voltages in terms of peak values, typically because of poor response time of voltage dividers. However the computed waveform shape in figure 7c is general as expected.

iii. The computed axial trajectory and speed, agree with experimental obtained time histories. Moreover, the behaviour with pressure, running the code at different pressures, agrees well with experimental results. The radial trajectories and speeds are difficult to measure. The computed trajectories figure 7e agrees with the scant experimental data available. The length of the radial structure is shown in figure 7f. Computed speeds radial shock front and piston speeds and speed of the elongation of the structure are shown in figure 7g.

iv. The computed inductance (figure 7h) shows a steady increase of inductance in the axial phase, followed by a sharp increase (rising by more than a factor of 2 in a radial phase time interval about 1/10 the duration of the axial phase for the NX2).

v. The inductive energy (0.5LI²) peaks at 70% of initial stored energy, and then drops to 30% during the radial phase, as the sharp drop of current more than offsets the effect of sharply increased inductance (figure 7i).

vi. In figure 7j is shown the work done by the magnetic piston, computed using force integrated over distance method. Also shown is the work dissipated by the dynamic resistance, computed using dynamic resistance power integrated over time. We see that the two quantities and profiles agree exactly. This validates the concept of half Ldot as a dynamic resistance.

vii. Dynamic resistance, DR (DR will be discussed in Module 2, Section 2.3; Note 2). The piston work deposited in the plasma increases steadily to some 12% at the end of the axial phase and then rises sharply to just below 30% in the radial phase. Dynamic resistance is shown in figure 7k. The values of the DR in the axial phase, together with the bank surge impedance, are the quantities that determine I_peak.

viii. The ion number density has a maximum value derived from shock-jump considerations, and an averaged uniform value derived from overall energy and mass balance considerations. The time profiles of these are shown in the figure 7l. The electron number density (figure 7m) has similar profiles to the ion density profile, but is modified by the effective charge numbers due to ionization stages reached by the ions.

ix. Plasma temperature too has a maximum value and an averaged uniform value derived in the same manner; are shown in figure 7n. Computed neon soft x-ray power profile is shown in figure 7o. The area of the curve is the soft x-ray yield in J. Pinch dimensions and lifetime may be estimated from figures 7e and 7f.

x. The model also computes the neutron yield, for operation in deuterium, using a phenomenological beam-target mechanism [25-27]. The model does not compute a time history of the neutron emission, only a yield number Y_n.

Thus as is demonstrated above, the model code when properly fitted is able to realistically model any plasma focus and act as a guide to diagnostics of plasma dynamics, trajectories, energy distribution and gross plasma properties.

1.4.4 Scaling parameters of the plasma focus pinch

The gross dynamics of the plasma focus is discussed in terms of phases. The dynamics of the axial and radial phases is computed using respectively a snowplow and an elongating slug model. A reflected shock phase follows, giving the maximum compression configuration of the plasma focus pinch. An expanded column phase is used to complete the post-focus electric current computation. Parameters of the gross focus pinch obtained from the computation, supplemented by experiments may also be summarised as follows:

Table 3

Plasma Focus Pinch Parameters	Deuterium	Neon (for SXR)
minimum radius rmin	0.15a	0.05a
max length (hollow anode) z	1.5a	1.6a
radial shock transit tcomp	5x10^-6a	4x10^-6a
pinch lifetime tp	10^-6a	10^-6a

where, for the times in s, the value of anode radius, a, is in m. For the neon calculations radiative terms are included; and the stronger compression (smaller radius) is due to thermodynamic effects.

1.5 Insights on plasma focus from numerical experiments using Lee model code

Using the Lee model code, series of experiments have been systematically carried out to look for behaviour patterns of the plasma focus.

Insights uncovered by the series of experiments include:

i. pinch current limitation effect as static inductance is reduced;

ii. neutron and SXR scaling laws;

iii. a global scaling law for neutrons versus storage energy combining experimental and numerical experimental data; and

iv. the nature and a fundamental cause of neutron saturation.

These significant achievements are accomplished within a period of twenty months of intensive numerical experimentation in 2008 to 2009. The numerical experimental research continues in 2010 with widening international collaboration.

References

[1] Lee, S. Plasma focus model yielding trajectory and structure. In Radiations in Plasmas, McNamara, B., Ed.; World Scientific, Singapore: 1984; Volume II, pp. 978–987.

[2] Glenn Millam Focus Fusion Society

[3] Muhammad Shahid Rafique. PhD thesis (in preparation) NTU, Singapore.

[4] S.Lee, S.P.Moo, C.S.Wong, A.C.Chew. Twelve Years of UNU/ICTP PFF- A Review. IC/98/231, ICTP Preprint, International Centre for Theoretical Physics, Trieste, Italy (1998), 101 pages.

[5] Lee S Radiative Dense Plasma Focus Computation Package: RADPF http://www.intimal.edu.my/school/fas/UFLF/File1RADPF.htm http://www.plasmafocus.net/IPFS/modelpackage/File1RADPF.htm

[6] S.Lee, T.Y.Tou, S.P.Moo, M.A.Eissa, A.V.Gholap, K.H.Kwek, S.Mulyodrono, A.J. Smith, Suryadi, W.Usala, M. Zakaullah. A simple facility for the teaching of plasma dynamics and plasma nuclear fusion. Amer J. Phys., USA, 1988, 56: 62-68.

[7] http://www.intimal.edu.my/school/fas/UFLF/

[8] S Lee and A Serban A, “Dimensions and lifetime of the plasma focus pinch,” IEEE Trans. Plasma Sci., 24, no.3, 1996, 1101-1105.

[9] Lee S. Diagnostics and Insights from Current waveform and Modelling of Plasma Focus. Keynote address: IWPDA, Singapore 2-July 2009

[10] Potter, D. E. The formation of high-density z-pinches. Nucl. Fusion. 1978, 18, 813–823.

[11] M.H.Liu. Soft X-Rays from Compact Plasma Focus. PhD thesis. NTU, Singapore, 1997

[12] S.Lee, P.Lee, G.Zhang, X. Feng, V.Gribkov, M.Liu, A.Serban & T.K.S. Wong. High Repetition High Performance Plasma Focus as a Powerful Radiation Source- IEEE Trans Plasma Science 26(4) 1119-1126 (1998).

[13] S. Lee, 2000–2007. [Online]. Available: http://ckplee.myplace.nie.edu.sg/ plasmaphysics/

[14] S. Lee, ICTP Open Access Archive, 2005. [Online]. Available:

http://eprints.ictp.it/85/

[15] G.X.Zhang. Plasma Soft x-ray source for Microelectronics lithography. PhD thesis, NTU, Singapore 1999.

[16] S.Lee, V.Kudryashov, P.Lee, G.Zhang, A..Serban, X..Feng, M.Liu, and T.K.S.Wong. Lithography Using a Powerful Plasma Focus Soft X-ray Source. International Congress on Plasma Physics, Prague, Czech Republic, June 1998. To be published in Procs.

[17] D Wong, P Lee, T Zhang, A Patran, T L Tan, R S Rawat and S Lee, “An improved radiative plasma focus model calibrated for neon-filled NX2 using a tapered anode,” Plasma Sources Sci. Technol. 16, 2007, pp. 116-123.

[18] E P Bogolyubov, V D Bochkov, V A Veretennikov, L T Vekhoreva, V A Gribkov, A V Dubrovskii, Yu P Ivanov, A I Isakov, O N Krokhin, P Lee, S Lee, V Ya Nikulin, A Serban, P V Silin, X Feng and G X Zhang, “A powerful soft x-ray source for x-ray lithography based on plasma focusing” 1998 Phys. Scripta., vol. 57, 1998, pp. 488-494.

[19] V Siahpoush, M A Tafreshi, S Sobhanian and S Khorram, “Adaptation of Sing Lee’s model to the Filippov type plasma focus geometry,” Plasma Phys. Control. Fusion 47, 2005, pp. 1065-1072.

[20] S. Lee and S.H. Saw “Neutron scaling laws from numerical experiments,” J. Fusion Energy, 2008, 27, no. 4, pp. 292–295.

[21] S. Lee S. “Current and neutron scaling for megajoule plasma focus machines,” Plasma Phys. Control. Fusion, 2008, 50, no. 10, p. 105 005 (14pp).

[22] S. Lee and S.H. Saw “Pinch current limitation effect in plasma focus,” Appl. Phys. Lett., 2008, 92, no. 2, p. 021 503.

[23] S. Lee S, P. Lee, S H. Saw and R.S. Rawat, “Numerical experiments on plasma focus pinch current limitation,” Plasma Phys. Control. Fusion, 2008, 50, no. 6, 065 012 (8pp).

[24] V. A. Gribkov, A. Banaszak, S. Bienkowska, A.V. Dubrovsky, I. Ivanova-Stanik, L. Jakubowski, L. Karpinski, R.A. Miklaszewski, M. Paduch, M.J. Sadowski, M. Scholz, A. Szydlowski and K. Tomaszewski “Plasma dynamics in the PF-1000 device under fullscale energy storage: II. Fast electron and ion characteristics versus neutron emission parameters and gun optimization perspectives,” J. Phys. D,Appl. Phys., 2007, 40, no. 12, 3592–3607

[25] S Lee, S H Saw, P C K Lee, R S Rawat and H Schmidt, “Computing plasma focus pinch current from total current measurement,” Appl. Phys. Lett. 92 , 2008, 111501

[26] Akel M., Al-Hawat Sh., Lee S. “Pinch Current and Soft x-ray yield limitation by numerical experiments on Nitrogen Plasma Focus”. J Fusion Energy DOI 10.1007/s10894-009-9238-6. First online 21 August 2009

[27] Saw S. H., Lee P. C. K., Rawat R. S. & Lee S. 2009 ‘Optimizing UNU/ICTP PFF Plasma Focus for Neon Soft X-ray Operation’ IEEE Trans on Plasma Sc, 2009, 37, 1276-1282.

[28] Saw S. H. and Lee S. ^“ Scaling laws for plasma focus machines from numerical experiments”. Invited paper: IWPDA, Singapore 2&3 July 2009

[29] Saw S. H. and Lee S. “Scaling the plasma focus for fusion energy considerations”. Tubav Conferences: Nuclear & Renewable Energy Sources, Ankara, Turkey, 28 & 29 September 2009.

[30] Lee S.”Nuclear fusion and the Plasma Focus”, Invited paper Tubav Conferences: Nuclear & Renewable Energy Sources Ankara, Turkey, 28 & 29 September 2009.

[31] Lee S., Saw S. H., Lee P. & Rawat R. S., “Numerical Experiments on Neon plasma focus soft x-rays scaling”, Plasma Physics and Controlled Fusion, 2009, 51, 105013 (8pp).

[32] Akel M., Al-Hawat Sh., Lee S. “Numerical Experiments on Soft X-ray Emission Optimization of Nitrogen Plasma in 3 kJ Plasma Focus SY-1 Using Modified Lee Model”, J Fusion Energy DOI 10.1007/s10894-009-9203-4 First online, May 19, 2009.

[33] Lee S., Rawat R. S., Lee P., S H Saw S. H., Soft x-ray yield from NX2 plasma focus, JOURNAL OF APPLIED PHYSICS, 2009, 106, 023309.

[34] Lee S. “Neutron Yield Saturation in Plasma Focus-A fundamental cause”

APPLIED PHYSICS LETTERS, 2009, 95, 151503 published online 15 October 2009

[35] Lee S., Saw S. H., Soto L., Moo S. P., Springham S. V., Numerical experiments on plasma focus neutron yield versus pressure compared with laboratory experiments, Plasma Phys. Control. Fusion, 2009, 51 075006

[36] Lee, S (2004) Characterising the Plasma Focus Pinch and Speed Enhancing the Neutron Yield. In: First Cairo Conference on Plasma Physics & Applications. International Cooperation Bilateral Seminars (Vol 34). Forschungszentrum Juelich GmbH, Juelich, Germany, pp. 27-33. ISBN 3-89336-374-2 http://eprints.ictp.it/39/

Part 1 Basic Course

Module 2: Universal Plasma Focus Laboratory-The Lee model code

Follow the instructions (adapted to EXCEL 2003) in the following notes. You may also wish to refer to the supplementary notes SP1.doc. Instructions are given in some details in order to accommodate participants who may not be familiar with EXCEL. Those who find the instructions unnecessarily detailed may wish to skip the unnecessary lines. The code seems to run unreasonably sometimes agonizingly slow when used with EXCEL 2007. So EXCEL 2003 is strongly preferred.

Summary

2.1 Introduction to the Worksheet

2.2 Configuring the Universal Plasma Focus Laboratory (UPFL)

2.3 Firing a shot in NX2

2.4 Studying the results

2.5 Exercise 1: Interpreting and recording data from the Worksheet

2.6 Conclusion

The material

You should have RADPFV5.15de.xls (contained in the e-folder “Code and Data” accompanying this file) on your Desktop for the next step. Please also ensure you have kept an identical original copy in a RESERVE folder. You are going to work with the desktop copy; and may be altering it. Each time you need an unaltered copy; you may copy from the reserve folder and paste it onto the desktop.

2.1 Introduction to the Worksheet

2.1.1 Opening the worksheet

(Note: Click means the ordinary click on the left button of the mouse; as distinct from the term Right Click, which means the special click on the right button of the mouse.)

Double click on RADPFV5.15de.xls

Work sheet appears and should look like Figure 1 [shown only as an example]; for the following please refer not to Figure 1 but to your worksheet.

Security pop-up screen appears.

Click on enable macros

[2007 Security Warning “Macros have been disabled” appears at top left hand corner of Worksheet with side box “options”.

Click on “options”- select the button “Enable this content”- click OK

After this procedure, the worksheet is macro-enabled and is ready for firing. ]

Figure 1. Appearance of worksheet-EXCEL 2003 version; EXCEL 2007 version should not look too different.

2.1.2 Preliminary orientation for setting controls

(For the following instructions, use your Excel Sheet; not the above image)

Device configuration:

(Note: Each Cell of the Excel Worksheet is defined by a Column alphabet A, B, or C….. and a Row number 1, 2 or 3 etc. The Column alphabets are shown along the top border of the worksheet. The Row numbers are shown along the left border of the Worksheet. For example, Cell A4 is located at column A row 4. Another example: A4-F9 refers to the block of cells within the rectangle bordered by row A4-F4, column A4-A9, row A9-F9 and column F4-F9; the larger orange-red bordered rectangle, containing 6x6 cells, near the top left of figure1.)

Locate Cells A4 to F9. These cells are for setting bank parameters, tube parameters, operating parameters and model parameters.

Taper: Control Cells for anode taper are normally inactivated by typing 0 (number zero) in Cell H7. Ensure that H7 is filled with 0 (number zero); unless anode taper feature is needed.

One Click Device: This control cell R4 allows choice of a specific plasma focus using numbers; currently 3 machines are available chosen with numbers 1, 2 or 3. Ensure that R4 is filled with the number 0. (Otherwise the code will keep defaulting to the selected machine ‘1’ or ‘2’ or ‘3’.)

2.1.3 Preliminary orientation for computed results

Cells A10-G13: computed characteristic quantities of the configured plasma focus.

Cells K6-M7: computed neutron yield, component & total; if operated in deuterium

Cell N6-N7: computed SXR line radiation

Cells H10-N11: computed durations of axial phase, radial phase and pinch phase and end time of radial phase.

Cells A15-AI17: dataline: contains data on row 17 with corresponding labels (and units) in rows 15 and 16. Data: E₀, RESF, c=b/a, L₀, C₀, r₀, b, a, z₀, V₀, P₀, I_peak, I_pinchstart, T_pinchmin, T_pinchmax, peak v_a, peak v_s, peak v_p, a_min (which is r_min), z_max, pinch duration, V_max, n_ipinchmax, Y_n, Qsxr, Qsxr%, f_m, f_c, f_mr, f_cr, EINP, taxialend, SF, ID and Q_line; others may be added from time to time.

This is a recently introduced very useful feature; enables computed data for each shot to be copied and pasted onto another sheet; so different shots may be placed in sequence, and comparative charts may be made.

Columns A20 to AP20: computed point by point results (data are correspondingly labeled in row A18 with units in row 19) for the following quantities respectively:

Time in ms, total current, tube voltage, axial position, axial speed, time of radial phase in ms measured from the start of axial phase, time of radial phase in ns from the start of the radial phase, corresponding quantities of current, voltage, radial shock position, radial piston position, radial pinch length, radial shock, piston and pinch elongation speeds, reflected shock position, plasma temperature, Joule power, Bremsstrahlung, recombination, line emission powers, total radiation power, total power, Joule, Bremsstrahlung, recombintion, line emission energies, total radiation energy, total energy, plasma self-absorption correction factor, black-body power, specific heat ratio and effective charge number, number thermonuclear neutrons, number beam target neutrons, number total neutrons, ion density, volume radiation power, surface radiation power, plasma self-absorption correction factor , radial phase piston work in % of E₀, neon SXR energy emission.

Each computed quantity as a function of time (displayed in the relevant column) is displayed in a column. After a run each of these columns is typically filled to several thousand cells.

Computed results are also summarized in 8 figures:

Figure 1: (Top left) total discharge current and tube voltage

Figure 2: (Top right) axial trajectory and speed

Figure 3: radial trajectories

Figure 4: total tube voltage during radial phase

Figure 5: radial speeds

Figure 6: plasma temperature

Figure 7: Joule heat and radiation energies

Figure 8: Joule power and radiation powers

An additional figure 8a on the right displays the specific heat ratio and effective charge number during the radial phase.

2.2 Configuring the Universal Plasma Focus Laboratory (UPFL)

2.2.1 Configuring the worksheet for a specific machine

As a first exercise we configure the UPFL so it operates as the NX2, the High-repetition rate neon focus developed for SXR lithography in Singapore.

The parameters are:

Bank: L₀=20 nH, C₀=28 mF, r₀=2.3 mW

Tube: b=4.1 cm, a=1.9 cm, z₀=5 cm

Operation: V₀=11 kV, P₀=2.63 Torr, MW=20, A=10, At-Mol=1 (these last 3 defines neon for the code i.e. molecular (atomic) weight, atomic number and whether atomic or molecular)

Model: massf (f_m)=0.0635, currf (f_c)=0.7, massfr (f_mr) =0.16, currfr (f_cr)=0.7; these are the mass and current factors for the axial and radial phases.

(Note: 1. These model parameters had been fitted earlier by us so that the computed total current best fits a measured total current trace from the NX2.

Note:2. We will carry out exercises in fitting model parameters in Module 3)

Configuring: Key in the following: (e.g. in Cell A5 key in 20 [for 20nH], Cell B5 key in 28 [for 28mF] etc.

A5 B5 C5 D5 E5 F5

20 28 4.1 1.9 5 2.3

Then A9 B9 C9 D9 E9

11 2.63 20 10 1

Then A7 B7 C7 D7

0.0635 0.7 0.16 0.7

You may of course find it easier to follow the labels in A4-F4, to key in A5-F5 for the relevant parameters; i.e. A4 states L₀ nH; so fill in below it in A5 20 ; and so on.

For identification purposes key in at B3 ‘NX2’

2.3 Firing a shot in NX2

Place the cursor in any blank non-active space, e.g. G8. (point the cursor at G8 and click the mouse). Press ‘Ctrl’ and ‘A’ (equivalent to firing a shot).

The programme runs and in less than a minute [provided you are running on EXCEL 2003, For some reason for EXCEL 2007, the process may take longer, sometimes much much longer] the run has completed and your worksheet will look something like figure 2 below:

Figure 2. Appearance of Worksheet after a shot.

Figure 3. Plasma focus current is distorted from unloaded current waveform.

In figure 3 is superimposed a current waveform (in blue; you do not have this waveform) of the plasma focus short-circuited across its input end insulator; with the current waveform (pink) you have just computed [see your worksheet figure 2] (In a later session you will learn how to do the short-circuit computation and superimposition).

Notes:

Note 1

The first important point to stress (and one that should never be forgotten) is that the plasma focus current waveform is very much distorted from the damped sinusoid of the L-C-R discharge without the plasma focus load (figure 3). The ‘distortions’ are due to the electrodynamical effects of the plasma motion, including the axial and radial dynamics and the emission of SXR from the Neon plasma. The way we use the code is based on the premise that the features of these ‘distortions’ contain the information of the plasma electrodynamics.

The plasma focus loads the electrical circuit in the same manner as an electric motor loads its driving circuit. The loading may be expressed as a resistance. More specifically we may compute the loading or ‘dynamic’ resistance as follows in Note 2; which shows that the dynamic resistance due to the motion in the axial phase is more than the stray resistance of the capacitor bank in the case of the NX2. Note 3 shows further that the dynamic resistance due to the plasma motion in the radial phase is so large as to completely dominate the situation. This causes the large current dip as shown in figure 3.

Note 2

As an example we may estimate the effect of one of the electrodynamical effects. The quantity (1/2)(dL/dt) is a dynamic resistance.

In the axial phase L=(m/2p)*ln(b/a)*z where m is permeability and z is the position of the current sheath.

Differentiating, 0.5*dL/dt= 10^-7 *ln(4.1/1.9)*axial speed~0.8 mW per 10⁴ m/s axial speed; or 0.8 mW per unit speed of cm/ms. At the peak axial speed of 6.6 cm/ms (see figure 2 of worksheet), that gives us a circuit loading of ~ 5 mW; which is reduced to 3.5 mW when we consider the effect of the current factor. This is more than the loading of the stray resistance r_o of 2.3 mW. So the axial motion of the current sheath is an important loading to the circuit.

Note 3

Continuing along this vein we may estimate the dynamic resistive loading of the current sheath motion in the radial phase when L = (m/2p)*ln(b/r_p)*z_f, where r_p= radial piston position and z_f = length of the elongating column; both r_p & z_f changing with time.

Thus dL/dt = (m/2p)*ln(b/r_p)*dz_f/dt +(m/2p)*z_f*(dr_p/dt)/r_p

= 2*10^-7*(ln(b/r_p)*dz_f/dt +z_f*(dr_p/dt)/r_p) [both terms RHS are positive]

In the section (2.4) below we will get from the output figures of the worksheet the following values at around the time of peak piston speed:

r_p~2.4 mm, z_f~15 mm, dr_p/dt~13.5 cm/ms [1.35*10⁵ m/s]; dz_f/dt~1.7*10⁵ m/s; Substituting into expression above, we get at the time of peak piston speed dL/dt~190 mW ; giving us (after considering current factor of 0.7) still around 130 mW of dynamic resistive loading due to the current sheath motion. This dynamic resistance (compared to r₀ of just 2.3 mW) dominates the current profile at this stage.

Note 4

d[LI]/dt generates an induced voltage; with one important component in this situation being I*(dL/dt). Since we have already estimated that dL/dt~0.19 W; multiplying this by 0.7x200 kA of current (which is the approx value of current at this time) gives us just under 30 kV. So we note that the dynamics at this time (just as the radial shock is going on axis) contributes a back voltage of ~30 kV through this term.

The other term L*(dI/dt) terms is negative; so the maximum induced voltage is considerably less than 30 kV, as you can see from figure 2.

Note 5

As a separate exercise which you may like to do one day: What is the basis for saying that (1/2)*(dL/dt) is a dynamic resistance? Can you show this by examining the power term in the situation when an inductance is changing? Compare the inductive power flow: (d/dt)(0.5*L*I²) and the total power flow: VI=I*(d/dt)(L*I).

What do you notice?

2.4 Studying the results

(The results are obtained from your Excel Sheet; not from the above images in figure 2)

Remember we are operating a neon plasma focus.

Here are some important quantities obtained from the data line in row 17.

Computed I_peak: L17 322 kA

I_pinch: M17 162 kA (pinch current at start of

pinch phase)

Peak tube voltage: V17 26.1 kV

k_min: S17 0.075 (r_min or a_min/a)

[you may also check this against figure 3 of the worksheet.]

Durations: H11-N11

Axial phase ends at 1.172 ms

Radial phase ends at 1.407 ms (add 1.172 to 0.235 ms) of which the last 26.2 ns is the pinch phase.

Now we study the various figures displayed on the worksheet, Sheet1 (also shown in figure 2).

Fig 1

Computed current trace; One point of interest is to locate the ends of axial and radial phases on this trace; as well as the start and end of the pinch phase. To do this, select Fig 1 (by pointing cursor on figure 1 and clicking). Then point cursor arrow at trace near peak and move until point 1.17 ms appears; that is the end of axial phase which is also the start of the radial phase.

Note: This point occurs not at the apparent start of current dip, but a little before that. There is no distinct indication on the trace that precisely marks this point. The term rollover may be a better term suggesting a smooth merging of the axial and radial phase. The apparent current dip occurs a little after the end of the axial phase.

Next locate point 1.41 ms which is the end of the radial phase. Also locate the point 1.38 ms which is the start of the pinch phase. There is no clear indication on the trace to mark this point either.

Fig 2

Select Fig 2 (with cursor) and read off the pink curve that the peak axial speed reached is 6.6 cm/ms. Confirm this on the data line; Cell P17. How many km per hour is this? And what is the Mach Number? 1 cm/ms=36,000 km/hr; so 6.6 cm/ms=237,600 km/hr

Expressing this speed in km/hr is to give an idea of how fast the speeds are in the plasma focus; it should give the idea also of temperature, since for strong shock waves (high Mach number motion) there is efficient conversion of energy from directed to thermal, i.e. from high kinetic energy to high temperature.

Mach number=speed/sound speed; finding this number is one of the questions for Exercise 1 (see below).

Fig 3

Select Fig 3. Read from dark blue curve that piston hits axis (radius=0) at 178ns after start of radial phase; and outgoing reflected shock (light blue) hits incoming piston (pink curve) at 210ns at radius of 2.1mm. The pinch phase starts at this 209ns and ends at 235 ns at a further compressed radius of 1.42mm.

Note the square of the ratio of pinching a/r_min is a measure of the how much the ambient density has been increased by the pinching effect.

Fig 4

Computed waveform of tube voltage during radial phase. Note the peak value of the tube voltage induced by the rapid plasma motion.

Fig 5

Select Fig 5. Note from the dark blue curve that peak radial shock speed is 20.4 cm/ms just before the radial shock hits the axis at 178 ns after start of radial phase. Also read from the pink curve that peak piston speed is 14.2 cm/ms reached just before the radial shock reaches its peak speed. Yellow curve shows column elongation speed. Note that these peak speeds are also recorded in the data line.

Other figures:

Select Fig 6: and read the peak temperature reached.

Select Fig 7: and read the various energies.

Select Fig 8: and read the various powers

Note that more charts are plotted on Sheet2 of RADPFV5.15de.xls These charts form a more complete picture of the plasma focus pinch, and may be used as starting guides for laboratory measurements of the various plasma properties.

2.5 Exercise 1: Interpreting and recording data from the worksheet.

Fill in the following blanks:

Q0: Given the speed of sound in neon at room temperature is 450 m/s (1600 km/hr), the Mach number of the peak axial phase speed, and of the peak radial phase speed (radially inwards shock speed) are _______ and _______. (Note: The peak axial speed can be found from Fig.2, and the peak radial speed can be found from Fig.5, for this particular plasma focus operation, also recorded in the dataline)

Q1: The peak temperature reached is _______K.

Q2: At that temperature the effective charge number (from small figure) is _______ and specific heat ratio has a range as follows _______.

Q3: There is a moment in time when the temperature jumps by a factor of approximately 2. This is at _______ ns from start of radial phase (Note: this happens at reflected shock according to the model)

Q4: Joule heating reached a maximum value of_______ J.

Q5: Total radiation reached a maximum value of_______ J (Note: the – sign for the radiation energy indicates energy taken out of the plasma by emission; ignore the – sign for this measurement)

Q6: Line Radiation reached a maximum value of _______ J.

Q7: Peak radiation power reaches a value of _______ W.

2.6 Conclusion

We had an introduction to the Worksheet of RADPFV5.15de.xls

We configured the UPFL as the NX2 at 11 kV 2.6 Torr neon.

We used properly fitted model parameters. (Note: Fitting model parameters will be covered in a future session).

We noted that the current waveform is distorted from damped sinusoid-like waveform (damped sinusoid-like waveform is the current waveform when the plasma focus is short-circuited).

We studied the computed results, including total current, tube voltage, pinch current, radial and axial trajectories, radial and axial speeds, plasma temperature and plasma Joule heating and radiation energies.

We also located various points on the current trace including: end of axial phase/start of radial phase; end of radial phase; start and end of pinch phase.

Note: This particular numerical ‘shot’ used properly fitted model parameters. The results of dynamics, electrodynamics and radiation as seen above are, in our experience, comparable with the actual experiments conducted at NTU/NIE.

End of Module 2.

Reference to the Lee model code should be given as follows:

Lee S. Radiative Dense Plasma Focus Computation Package (2011): RADPF www.plasmafocus.net http://www.intimal.edu.my/school/fas/UFLF/

Part 1 Basic Course

Module 3: I. Configuring and fitting computed current to measured current

II. Comparing a large PF with a small PF-neutron yield etc

(Follow the instructions in the following notes. You may also wish to refer to the supplementary notes SP2.doc)

Summary

For this module we fit model parameters so that computed current waveform matches measured current waveform.

First we configure the UPFL (RADPFV5.15de.xls) for PF1000 operating in deuterium; using trial model parameters. We fire a shot. We do not know how good our results are without a reference point; i.e. some comparison with experimental results.

A total current waveform of the PF1000 has been published; we have it in digitized form in a file PF1000data.xls. We also have the chart of this waveform displayed in this file.

To ensure that our computed results are comparable to experimental results, the key step is to fit model parameters, by adjusting the model parameters until the computed total current trace matches the measured total current trace.

To do this, we add PF1000data.xls to our numerical focus laboratory RADPFV5.15de.xls. Next we plot the computed current waveform in the same chart. The model parameters are varied; at each variation the focus is fired, and the computed current waveform is compared with the measured waveform. The process is continued until the waveforms are best matched. A good match gives confidence that the computed results (trajectories, speeds, temperature, neutron and radiation yields etc) are comparable with actual experimental results.

After the guided fitting of the PF1000, we have a self exercise to fit the Chilean PF400J.

We then tabulate important results of both machines, and do a side-by-side comparison of big versus a small plasma focus to obtain important insights into scaling laws/rules of the plasma focus family.

Configuring and fitting computed current to measured current

3.1 Configure the code for the PF1000 using trial model parameters

3.2 Place a measured (published) PF1000 current waveform on Sheet3

3.3 Place the computed current waveform onto the same chart as the measured current waveform in Sheet3

3.4 Vary the model parameters to obtain matching of computed versus measured current traces

II. Comparing a large PF with a small PF-neutron yield etc

3.5 Exercise 2: Tabulate results for PF1000 obtained in numerical experiments

3.6 Exercise 3: Fitting the PF400J and tabulate the results for PF400J side by side with the results for PF1000, for a comparative study

3.7 Conclusion

The material

You need RADPFV5.15de.xls for the following work. Copy and paste a copy on your Desktop. You also need the files PF1000data.xls, PF400data.xls and compareblank.xls for this session. Copy and paste a copy of each file onto your desktop.

I. Configuring and fitting computed current to measured current (guided–4 hrs)

3.1 Configure the code for the PF1000 using trail model parameters

Double click on RADPFV5.15de.xls on your Desktop.

Security pop-up screen appears.

Click on enable macros

The worksheet opens

[EXCEL 2007: Security Warning “Macros have been disabled” appear at top left hand corner of Worksheet with side box “options”.

Click on “options”- select the button “Enable this content”- click OK]

After this procedure, the worksheet is macro-enabled and is ready for firing.

Type in cell B3: PF1000; for identification purposes.

The PF1000, at 40 kV, 1.2 MJ full capacity, is one of the biggest plasma focus in the world. Its 288 capacitors have a weight exceeding 30 tonne occupying a huge hall. It is the flagship machine of the ICDMP, International Centre for Dense Magnetised Plasmas.

We are now working on this Plasma Focus.

We use the following bank, tube and operating parameters for the PF1000:

Bank: L₀=33.5 nH, C₀=1332 mF, r₀=6.3 mW

Tube: b=16 cm, a=11.55 cm, z₀=60 cm

Operation: V₀=27 kV, P₀=3.5 Torr, MW=4, A=1, At-Mol=2

For this exercise we do not know the model parameters. We will use the trial model parameters recommended in the code (See cells T9-V9)

Model: massf (f_m) = 0.073, currf (f_c )= 0.7, massfr (f_mr) = 0.16, currfr (f_cr) = 0.7; our first try.

Configuring: Key in the following:

A5 B5 C5 D5 E5 F5

33.5 1332 16 11.55 60 6.3

Then A9 B9 C9 D9 E9

27 3.5 4 1 2

Then A7 B7 C7 D7

0.073 0.7 0.16 0.7 for first try

Or follow the guide in A4-F4, to key in A5-F5 for the relevant parameters.

Fire a shot: Place the cursor in any blank non-active space, e.g. G8. (point the cursor at G8 and click the mouse). Press ‘Ctrl’ and ‘A’. (firing a shot)

The program runs and results are displayed in columns and also in figures.

Is our simulation any good? Not if there is no reference point!!

To assess how good our simulation is, we need to compare our computed current trace with the measured current trace, which has been published.

Note that at this point: the configured RADPFV5.15de.xls contains computed data for PF1000 with the trial model parameters of: massf (f_m) = 0.073, currf (f_c) = 0.7, massfr (f_mr) = 0.16, currfr (f_cr) = 0.7

3.2 Place a measured (published) PF1000 current waveform on Sheet3

3.2.1 The PF1000 current waveform is in the file PF1000data.xls. You now want to place this data file as an additional sheet in our RADPFV5.15de.xls workbook.

RADPFV5.15de.xls is already open, click on ‘File’ tab; drop down appears,

click on ‘Open’.

Look in: Desktop; select PF1000data.xls; double click to open this file.

Click on ‘Edit’ tab; select ‘Move or Copy Sheet’.

A window pops out; ‘Move selected sheets To book’; select RADPFV5.15de.xls. ‘Before Sheet:’ select ‘(move to end)’.

Click ‘OK’.

You have copied PF1000data.xls into RADPFV5.15de.xls as Sheet1(2); you might like to rename Sheet1(2) as Sheet3.

[EXCEL 2007: ADPFV5.15de.xls is already open, minimize it by clicking on top right hand corner tab with the - sign.

Open PF1000data.xls.

Locate tab Sheet3 on lower left corner of worksheet.

Right click on tab Sheet3.

Select move or copy to book RADPFV5.15de.xls

Click on (move to end)

Tick Create a copy

Click OK]

3.2.2 With this procedure you have copied PF1000data.xls as Sheet3 in RADPFV5.15de.xls. The chart has already been prepared. The measured current waveform appears in the chart.

3.3 Place the computed current waveform onto the same chart as the measured current waveform in Sheet3

In the next steps we will place the computed current data from Sheet1 into this same chart in Sheet 3, by the following procedure.

Place the cursor on the chart; click, then right click; drop down appears- click on source data; click on series. In the series box click on computed current in kA; then in the box against ‘X Values’ type in the following string: “=sheet1!$a$20:$a$6000” [without the quotation marks]. Next click in the box against ‘Y Values’ and type in the following: “=sheet1!$b$20:$b$6000” [without the quotation marks]. Click button ‘OK’.

[EXCEL 2007: Position the cursor on the chart in Sheet3 containing the measured current waveform. Now right click. Pop-up appears. Click on Select data. Select the series “Computed current in kA”. Click on Edit. On the new pop-up for series the name “=”Computed current in kA” is already there.

For Series X values key in “=Sheet1!$a$20:$a$6000”

For Series Y values key in “=Sheet1!$b$20:$b$6000”

Click OK; and click OK.]

The computed current waveform from Sheet1 is charted in the figure in Sheet3 with the same time scale and the same current scale.

You can now compare the computed current trace with the measured current trace.

You should see a pink trace which has just appeared on the chart. The pink trace (see figure 1 below) is the computed current trace transferred from Sheet1 (where the time data in ms is in column A, from A20-A several thousand; and corresponding computed current data in kA is in column B, from B20 to B several thousands). We are selecting the first 5980 points (if that many points have been calculated) of the computed data; which should be adequate and suitable.

Figure 1. Comparison of traces: Note that there is very poor matching of the traces; using the first try model parameters.

3.4 Vary the model parameters to obtain matching of computed versus measured current traces

Note that bank, tube and operating parameters have all been given correctly.

3.4.1 First fit the axial phase

[suggestion: read SP2.doc pg 3 ‘First step is fitting the axial phase’.]

From the comparison chart on Sheet2,

We note: that the computed current dip comes much too early;

that the computed current rise slope is only very slightly low;

that the computed current maximum is too low.

All these 3 observations are consistent with a possibility that the axial speed is too fast; which would cause the radial phase to start too early. Too high an axial speed would also cause too much loading on the electrical circuit (similar to the well known motor effect) as the quantity [0.5 x dL/dt=0.5x L’x dz/dt] is a dynamic resistance loading the circuit during the axial phase; here the inductance per unit length L’=(m/2p) x ln(b/a). This too high speed would also lower the peak current.

To reduce the axial speed, we could increase the axial mass factor. We note that the axial phase ends too early by some 20%; indicating the axial speed is too fast by 20%.

In the plasma focus (as in pinches, shocks tubes and other electromagnetically driven plasma devices) speed~density^0.5. So the correction we need is to increase the axial mass factor by 40%. So try an axial mass factor of 0.073x1.4~ 0.1.

We toggle to Sheet1 by clicking on ‘Sheet1’ (just below the worksheet).

Click on cell A7, and type in 0.1.

Fire the focus by pressing Ctrl+A.

Program runs until completed, and results are presented.

Note TRadialStart (H11) has increased some 0.6 ms.

Toggle to Sheet3 (i.e. click on Sheet3 just below work sheet).

Note that the computed current dip is now closer to the measured current dip in time (still short by some 10%; reason being that increasing the axial mass factor reduces the speed which in turn causes a reduced loading. This increases the current which tends to increase the axial speed so that our mass compensation of 40% becomes insufficient). The value of the computed peak is also closer to the measured. So we are moving in the right direction!

But still need to move more in the same direction. Next try axial mass factor of 0.12. Toggle to Sheet1, type 0.12 in A7. Fire. Back to Sheet3. Note improvement in all 3 features.

In similar fashion, gradually increase the axial mass factor. When you reach 0.14 you will notice that the computed current rise slope, the topping profile, the peak current and the top profile are all in good agreement with the measured. The computed trace agrees with the measured trace up to the start of the dip. Note that the axial model parameters at this stage of agreement are: 0.14 and 0.7. You may wish to try to improve further by making small adjustments to these parameters. Or else go on to fit the radial model parameters.

3.4.2 Next, fit Radial phase

Note that the computed current dip is too steep, and dips to too low a value. This suggests the computed radial phase has too high a speed. Try increasing the radial mass factor (cell C7), say to 0.2. Observe the improvement (dip slope becomes less steep) as the computed current dip moves towards the measured dip.

Continue making increments to mass(f_r) (cell C7). When you have reached the mass(f_r) value of 0.4; it is becoming obvious that further increase will not improve the matching; the computed dip slope has already gone from too steep to too shallow, whilst the depth of the dip is still excessive.

How to raise the bottom of the dip? Here we suppose the following scenario:

Imagine if very little of the current flows through the pinch, then most of the total current will flow unaffected by the pinch. And even if the pinch were a very severe one, the total current (which is what we are comparing here) would show hardly a dip. So reducing the radial current fraction, ie currfr (or f_cr) should reduce the size of the dip.

Let us try 0.68 in cell D7. Notice a reduction in the dip. By the time we go in this direction until currfr(f_cr) is 0.65, it becomes obvious that the dip slope is getting too shallow; and the computed dip comes too late.

One possibility is to decrease massfr(f_mr) (which we note from earlier will steepen the dip slope); which however will cause the dip to go lower; and it is already too low. Another possibility is to decrease the axial phase massf(f_m), as that will also move the computed trace in the correct direction.

Try a slight decrease in axial massf (f_m), say 0.13.

Note that this change aligns the dip better but the top portion of the waveform is now slightly low, because of the increased loading on the electrical circuit by the increase in axial speed. This suggests a slight decrease to circuit residual resistance r₀ (or changes to L₀ or C₀; fitting those could be tricky, and we try to avoid unless there are strong reasons to suspect these values). Easier to try lowering r₀ first. Try changing r₀ to 6.1 mW.

The fit is quite good now except the current dip could be steepened slightly and brought slightly earlier in time. Decrease massfr(f_mr), say to 0.35. The fit has improved, and is now quite good, except that the dip still goes too low. At this stage we check where we are at.

Toggle to Sheet1. Note from Sheet1 that the radial phase ends at 9.12 ms. Back to Sheet3.

Observe (using cursor) that the point 9.12 is not at the point where the computed (pink curve) dip reaches its inflection point; but some 0.02 ms before that point (see figure below).

So we note that the computed curve agrees with the measured curve up to the end of the radial phase with a difference of less than 0.02 MA out of a dip of 0.66 MA (or 3%).

The fitting has already achieved good agreement in all the features (slopes and magnitudes) of the computed and measured total current traces up to the end of the radial phase.

Do not be influenced by agreement, or disagreement of the traces beyond this end point.

Figure 2. The best fit.

So we have confidence that the gross features of the PF1000 including axial and radial trajectories, axial and radial speeds, gross dimensions, densities and plasma temperatures, and neutron yields up to end of radial phase may be compared well with measured values.

Moreover the code has been tested for neutron and SXR yields against a whole range of machines and once the computed total current curve is fitted to the measured total current curve, we have confidence that the neutron and SXR yields are also comparable with what would be actually measured.

Having said that, those of you who have some experience with the plasma focus would note that at the end of the radial phase, some very interesting effects occur leading to a highly turbulent situation with occurrence, for example, of high density hot spots. These effects are not as yet modeled in the code. Despite this drawback, the postulated beam-target neutron yield mechanism seems able to give estimates of neutron yield which broadly agree with the whole range of machines. For example, the neutron yield computed in this shot of 1.08x10¹¹ is in agreement with the reported PF1000 experiments.

One further note: We have recently confirmed that the above discussion of fitting applies typically to machines with low L₀, below perhaps 60 nH. For machines above 100 nH another strategy of fitting or even modelling may need to be adopted. This is related to the comment just above this note.

II. Comparing a large PF with a small PF - neutron yield etc

3.5 Exercise 2: Tabulate results for PF1000 obtained in numerical experiments

You have been following the guided steps in the above fitting:

Fill in the following:

Q1: My best fitted model parameters for PF1000, 27kV, 3.5 Torr deuterium are:

fm=

fc=

fmr=

fcr=

Q2: Insert an image of the discharge current comparison chart in Sheet3 here.

Q3: Fill up the following table. To help with collecting the data use the ‘dataline’ feature in Sheet1 which is the line 17 with corresponding descriptors in lines 15 & 16. To tabulate the following, use the file compareblank.xls already prepared for this purpose.

Parameter	PF 1000 (at 27 kV 3.5 Torr D2)
Stored Energy E0 in kJ
Pressure in Torr, P0
Anode radius a in cm
c=b/a
anode length z0 in cm
final pinch radius rmin in cm
pinch length zmax in cm
pinch duration in ns
rmin/a (rmin is also called amin)
zmax/a
Ipeak in kA
Ipeak/a in kA/cm
S=(Ipeak/a)/(P0^1/2)( kA/cm)/Torr^1/2
Ipinch in kA
Ipinch/Ipeak
Peak induced voltage in kV
peak axial speed in cm/ms
peak radial shock speed cm/ms
peak radial piston speed cm/ms
peak temperature in 10^6K
neutron yield in units of 10^6

[After filling, save this Excel sheet you will use the same Excel sheet to fill in the results for PF400J which is the subject of the next exercise.]

3.6 Exercise 3: Fitting the PF400J and tabulate results for PF400J side by side with the results for PF1000 for a comparative study

Participants are to fit computed current to measured current waveform of PF400J (bank, tube and operating parameters all correctly given)

In Module 2, we worked with the Singaporean NX2; a 3kJ neon plasma focus designed for SXR lithography. For our first fitting exercise we worked with the Polish PF1000, one of the largest plasma focus (MJ) in the world. You are now given data for the PF400J, a small sub-kJ plasma focus operated in Chile at the Atomic Energy Commission, for the specific purpose of investigating small focus devices. [Note: The 1000 in PF1000 refers to kJ. The 400 in PF400J refers to J; so the energy ratio of the PF1000 to the PF400J is: 1,200,000/400~3000]. The PF400J is a small table top device with all components fitting on a small table top. The PF1000 has a huge chamber and its capacitor bank fills a whole big hall.

Given: the current waveform data of the PF400J, digitized from a published waveform. The data is in the file PF400data.xls.

Your task: is to fit model parameters until the computed current waveform matches the measured waveform. Some guidance is given below.

Suggested steps to fit PF400J

i. Make a clean copy of RADPF05.15de.xls from your Reserve folder to your Desktop. Open this file.

ii. Copy PF400data.xls as Sheet3 of RADPF05.15de.xls using procedure as in section 3.2.1 above. The measured waveform is already pre-charted.

iii. Transfer computed current data from Sheet1 onto the measured current chart in Sheet3; as in step in section 3.3 above using strings: “=sheet1!$a$20:$a$6000” [without the quotation marks] and “=sheet1!$b$20:$b$6000” [without the quotation marks]. No trace of computed current appears yet, since we have not yet ‘fired’ PF400J.

Write down the bank, tube and operating parameters (from the table in the lower part of Sheet3, NOT from the top line, which contains some nominal values). Toggle to Sheet1.

iv. Configure the Universal Plasma Focus with the following bank, tube and operating parameters for the PF400

L0(nH)	C0(μF)	b(cm)	a(cm)	z0(cm)	r0(mW)
40	0.95	1.55	0.6	1.7	10
massf	currf	massfr	currfr	model parameters
V0(kV)	P0(Torr)	MW	At No.	At1;Mol2	Operation Parameters
28	6.6	4	1	2

Key in the first try model parameters; [scroll a little to the right and use the suggested parameters for the UNU ICTP PFF, cells T9-V9].

v. Fire PF400J; and see the comparative results by toggling to Sheet3.

vi. Fitting the computed current waveform to the measured waveform

vii. Suggested first steps: Fit the axial region by small adjustments to f_m and f_c, where necessary. In fitting the axial phase, the more important region to work on is the later part of the rising slope and the topping profile towards the end of the axial phase. So each time you should note the position of the end of the axial phase from Sheet1 and locate that position on the chart in Sheet3, using the cursor.

viii. Final steps: When you have done the best for the axial phase up to the end of the axial phase, then proceed to fit the radial phase. Tip: The dip for the PF400J is not very dramatic. Enlarge the trace so the rollover and the dip can be more clearly compared.

Fill in the following questions, copy and paste and e-mail to me.

Questions

Q1: My best fitted model parameters for PF400J, 28 kV, 6.6 Torr deuterium are:

fm=

fc=

fmr=

fcr=

Q2: Insert an image of the discharge current comparison chart in Sheet3 here.

Q3: Complete the Excel Sheet which you started in the last Exercise; to compare a BIG (~500 kJ) plasma focus with a small one (~400J). As you fill up, note particularly each group of ratios (each group is denoted by a different colour). Note particularly the order of magnitude of the ratios. [use the Excel sheet, rather than this table].

The ratios below were calculated from the actual PF1000 and PF400J results; and left here as a check for you. Calculate your own ratios from your own results. At the end of the exercise save this Excel Sheet as PFcomparison.xls. It will be used again if eventually you go to the more advanced exercises of Modules 5.

Make up the following table comparing a BIG plasma focus with a small plasma focus.

Parameter	PF1000 ( at 27kV 3.5 Torr D2)	Ratio PF1000/PF400J	PF400J (at 28kV 6.6 Torr D2)
Stored Energy E0 in kJ	486	1313	0.37
Pressure in Torr, P0	3.5	0.53	6.6
Anode radius a in cm	11.55	19.3	0.6
c=b/a	1.39	0.54	2.6
anode length zo in cm	60	35.2	1.7
final pinch radius rmin in cm		26.7
pinch length zmax in cm		22.2
pinch duration in ns		53
rmin/a		1.4
zmax/a		1.16
Ipeak in kA		14.6
Ipeak/a in kA/cm		0.76
S=(Ipeak/a)/(P0^1/2)( kA/cm)/Torr^1/2		1.05
Ipinch in kA		9.64
Ipinch/Ipeak		0.65
Peak induced voltage in kV		2.4
peak axial speed in cm/ms		1.24
peak radial shock speed cm/ms		0.48
peak radial piston speed cm/ms		0.48
peak temperature in 10^6K		0.19*
neutron yield Yn in 10^6		81920
Measured Yn in 10^6: range	(2 - 7)E+03		0.9-1.2
Measured Yn in 10^6 :highest	2.0E+04

We could then use the tabulation for several applications including the following:

Think of scaling rules, laws:

Q4: What is the significance of the Speed Factors S of PF1000 and PF400J? Which one’s temperature should be higher?

***The ratio radial speed/axial speed is:

Note: http://www.plasmafocus.net/; download the Theory of the model

3.7 Conclusion

In this module we have learned how to fit a computed current trace with a measured current waveform, given all bank, tube and operational parameters. For the PF1000 we obtained a good fit of all features from the start of the axial phase up to the end of the radial phases; giving confidence that all the computed results including trajectories and speeds, densities, temperatures and neutron and radiation yields are a fair simulation of the actual PF1000 experiment.

We also fitted the computed current to the measured current of the PF400J; thus computing its dynamics and plasma characteristics and neutron yield.

We tabulated important results of the two machines side by side.

We noted important physics:

Although the machines differ greatly in storage energy and hence in physical sizes, the speed factor S is practically the same. This has given rise to the now well-known observation that all plasma focus, big and small, all operate with essentially the same energy per unit mass when optimized for neutron yield. See e.g.: http://en.wikipedia.org/wiki/Dense_plasma_focus (The Wikipedia article was written by IPFS and added on by various other groups and researchers.)

The axial speed is also almost the same; in which case the radial speeds would have been almost the same, except they (the radial speeds) are influenced by a geometrical factor [(c²-1)/lnc]^0.5. For these 2 machines the factors differ by 1.5; hence explaining the higher radial speeds in PF400J; and also the higher temperatures in the smaller PF400J.

The pinch dimensions scale with ‘a’ the anode radius. The pinch duration also scales with ‘a’, modified by the higher T of the PF400J, which causes a higher small disturbance speed hence a smaller small disturbance transit time. In this model this transit time is used to limit the pinch duration.

Finally we may note that just by numerical experiments we are able to obtain extensive properties of two interesting plasma focus machines apparently so different from each other, one a huge machine filling a huge hall, the other a desk top device. Tabulation of the results reveals an all important characteristic of the plasma focus family. They have essentially the same energy per unit mass (S).

A final question arising from this constant energy/unit mass: Is this at once a strength as well as a weakness of the plasma focus?

End of Module 3.

Reference to the Lee model code should be given as follows:

Lee S. Radiative Dense Plasma Focus Computation Package (2011): RADPF www.plasmafocus.net http://www.intimal.edu.my/school/fas/UFLF/

Part 1 Basic Course

Module 4: PF1000 neutron yield versus pressure

(Follow the instructions in the following notes. You may also wish to refer to the supplementary notes SP3.doc)

Summary

This module looks at variation of neutron yield with pressure; running PF1000 from short- circuit (very high pressure), through optimum pressure to low pressure. The very high pressure of the short-circuit shot stops all current sheath motion thus simulating a short circuit. The aim of this shot is just to obtain short-circuit current waveform for comparison with the focusing waveforms at different lower pressures. When the operating (ambient) pressure is low enough neutrons are emitted. The variation of the yield and other properties with pressure are compiled together, and presented on one chart in normalized form. In this way correlation of various quantities may be seen.

4.1 Configure the code for the PF1000 at 27 kV, 3.5 Torr D₂ using model parameters which we had fitted earlier

4.2 Fire the PF1000 at very high pressure, effectively a short circuit

4.3 Fire the PF1000 at lower pressures from 19 Torr down to 1 Torr; looking for optimum neutron yield

4.4 Exercise 4: Place the current waveforms (from section 4.3) at different pressures on the same chart for comparative study

4.5 Exercise 5: Tabulate results at different pressures for comparative study; including speeds, pinch dimensions, duration, temperature and neutron yield

4.6 General notes on fitting, yield scaling and applications of the Lee model code

The material

You need RADPFV5.15de.xls for the following work. Copy and Paste on your Desktop. You also need the files PF1000pressureblank.xls. This file contains also tabulation blanks for your convenience.

4.1 Configure the code for PF1000 at 27 kV, 3.5 Torr deuterium using model parameters which we had fitted earlier

i. Preparing Sheet3

Open RADPFV5.15de.xls; copy PF1000pressureblank.xls as Sheet3; using procedure which we have practiced in Module 3.

Examine Sheet3

PF1000pressureblank.xls, copied as Sheet3 has the parameters of PF1000 recorded along the top rows. One set of measured time-current data is supplied at Columns A and B. To save participants some time, time-current data for several traces are already computed and filled in: 3.5 Torr at columns C and D; 19 Torr at columns G & H; 7.5 Torr at columns M & N and 1 Torr at columns S & T. These are already plotted in figure 1. You are required to fire shots at several other pressures and copy the time-current data for these shots and paste onto the columns reserved for these shots. In this way you fill up the charts with sufficient traces to cover the optimum pressures for neutron yield. These pressures are 100,000, 14, 10, 6 and 2 Torr.

Scrolling to the right you see table 1 with plasma focus properties at various pressures. The properties corresponding to the shots at 19, 7.5, 3.5 and 1 Torr are already filled in. Data from those other shots to be fired by the participants are to be filled in to complete the table. Below table 1 is table 2 with data normalized from table 1 using the shot at 7.5 Torr as the reference shot for normalizing.

Figure 2 displays the normalized Y_n, I_peak, I_pinch and radial phase piston work EINP versus pressures from table 2.

[Note that the curves in figure 2 have places where they all come to zero. At the start the curves lookvery messy and may appear strange. That is because the table has not been filled for those shots yet to be fired. For these shots the data points have been put to zero. The curves will take on their correct shapes once the data has been correctly filled in. This is the job for the participant]

ii. Configure the code for PF1000

Use the data in PF1000 pressureblank.xls to configure.

Bank: L₀=33.5 nH, C₀=1332 mF, r₀=6.1 mW

Tube: b=16 cm, a=11.55 cm, z₀=60 cm

Operation: V₀=27 kV, P₀= ? Torr, MW=4, A=1, At-Mol=2

Model: f_m=0.13, f_c=0.7, f_mr=0.35, f_cr=0.65

4.2 Fire the PF1000 at very high pressure, effectively a short circuit

Key in 100,000 Torr at B9.

[Note: In the laboratory it is of course impossible to fire such a shot and a physical short-circuit may need to be used at the insulator end of the plasma focus; or fire at the highest safe pressure in argon. In the lab we have used 50 Torr argon, to obtain very approximate results.]

[In the numerical experiment at this high pressure the current sheath only moves a little down the tube, adding hardly any inductance or dynamic loading to the circuit. So it is equivalent to short-circuiting the plasma focus at its input end. In the code there is a loop during the axial phase, computing step- by- step the variables as time is incremented. The loop is broken only when the end of the anode (non-dimensionalised z=1) is reached. In this case we do not reach the end of the anode. However there is an alternative stop placed in the loop that stops the run when time reaches nearly 1 full cycle is completed (non-dimensionalised time=6 ie nearly 1 full cycle time, 2p, of the short- circuited discharge).

At the start of the run, the code computes a quantity ALT= ratio of characteristic capacitor time to sum of characteristic axial & radial times. Numerical tests have shown that when this quantity is less than 0.65, the total transit time is so large (compared to the available current drive time) that the radial phase may not be efficiently completed. Moreover because of the large deviation from normal focus behaviour, the numerical scheme and ‘house keeping’ details incorporated into the code may become subjected to numerical instabilities leading to error messages. To avoid these problems a time-match guard feature has been incorporated to stop the code from being run when ALT<0.65. When this happens one can over-ride the stop; and continue running unless the run is then terminated by Excel for e.g. ‘over-flow’ problems. In that case one has to abandon the run and reset the code.]

Fire the high pressure shot. The pressure is too high for a normal run and we are automatically toggled over to the Macro; the Visual Basics Code appears at Statement 430 Stop; with a warning message that pressure is too high. In this case we know what we are doing, and over-ride as follows: Click on ‘Run’ (above the code sheet), and ‘continue’. Another ‘Stop’ appears just below line 485; with a warning about transit time. Click on ‘Run’ and ‘continue’; another ‘Stop’ appears below Line 488. Click on ‘Run’ and ‘Continue’.

In a little while the run has proceeded and finally the statement “If T > 6 Then Stop” appears; indicating we have completed nearly one cycle of the capacitor discharge; and reached the pre-set time limit.

Now, locate the ‘x’ at the extreme right hand corner of the screen. Click on this ‘x’; pop-up appears with the message ‘This command will stop the debugger’. Click on OK, which toggles us back to the worksheet, Sheet1.

Copy the data in columns A & B from A20 and B20 to the end of the computed current data (several thousand cells down); toggle to Sheet3 and ‘paste’ the copied time-current data onto Columns E & F (in the labeled space provided in Sheet3. Locate table 1 by scrolling to the right. Fill in the value of I_peak [read from figure 1 or from the relevant cell of the dataline] onto the table 1 against 100,000 Torr. Put zero against all the other quantities (I_pinch, peak v_a, S, peak v_s …. T and Y_n….)

4.3 Fire the PF1000 at lower pressures from 19 Torr down to 1 Torr, looking for optimum neutron yield

Fire the next shot at 14 Torr. As the ALT value is over 0.65, the run proceeds as normal. Copy the time-current data from Columns A & B (from rows 20 down) to Sheet3 columns I & J. Fill in the table 1 [I_peak, I_pinch, peak v_a, S, peak v_s … T… Y_n…n_i & EINP taken from the data line) for the data from shot 14 Torr.

Repeat for pressures 10, 9, 8, 7.5 , 7, 6, 3.5 , 2 and 1 Torr; tabulating the data for all these shots onto table 1; but copy and paste the time-current data for only selected shots of 14, 10, 6 and 2 [in order for figure 1 not to become too crowded]. The list of pressures had been chosen as above in order to demonstrate the way the neutron yield varies with pressure. It is clear that Y_n increases rapidly from 14 Torr to 10 Torr. More points are chosen between 10 Torr and 6 Torr and it is obvious that the optimum pressure (for Y_n) is between 8 Torr and 7 Torr. The participant will notice this as the shots are fired and as the Y_n data is copied on to the table 1.

4.4 Exercise 4: Place the current waveforms (from section 4.3) at different pressures on the same chart for comparative study

Suggested procedure: To save you time, the comparison chart, figure 1 has already been created for you, and pre-filled with several waveforms namely 19, 7.5, 3.5 and 1 Torr. You only have to fill in the ones for 100,000 and 14, 10, 6 and 2 Torr in the correct columns indicated by the column headings already placed on Sheet3.

You will note that the computed current waveform for 3.5 Torr falls neatly on the measured current waveform (as you have seen during an earlier exercise precisely with this PF1000 27 kV, 3.5 Torr current waveform.) You will recognize that we are using the computed 3.5 Torr shot to fit our UPFL to the PF1000 to obtain the model parameters for the PF1000.

Thereafter the assumption is that the model parameters apply for all the other shots.

It would of course be better if for every pressure or every shot we have a measured current trace to fit the code. However despite this assumption about the model parameters, the numerical experiment does show some very interesting features as we proceed below.

4.5 Exercise 5: Tabulate results at different pressures for comparative study; including speeds, pinch dimensions, duration, temperature and neutron yield

This tabulation has already been done as step (4.3) proceeded above.

In order to chart some of the computed data on one comparative chart, as mentioned already, as you fill in table 1, table 2 is at the same time filled in with each data column normalized to the data at 7.5 Torr, which was found to be the pressure with the highest Y_n. Thus the values of all the data in the normalized table is in the region of 1.

Plot normalized Y_n, I_peak, I_pinch, and radial EINP against P₀.

[As you fill in table 1, the normalized quantities are automatically computed, and the chart begins to take the correct shape. At the start the chart is in a jumble because many points have not been filled in, and thus there are erratic zero points all over the place.]

Discussion

Note 1

Look at the change of current waveforms from very high pressures to low pressures. At very high pressures the waveform is a damped sinusoid. At 19 Torr the characteristic flattening of the current waveform due to dynamics is already clearly evident. The current peak comes earlier and is lower than the unloaded (high pressure) case, the current then droops until the rollover into the dip (due to the increased radial phase loading) at around 15 μs. At lower pressures these characteristics remain the same except that the current trace is depressed more and more as speed increases. The peaking (reaching maximum current) also comes earlier and earlier, as does the radial phase rollover of the current trace. This is characteristic of an L-C-R circuit with increasing resistance R, as R increases from light towards critical damping.

At 2.6 Torr, there is hardly any droop, the current waveform showing a distinct flat top leading to the rollover. At 1 Torr the axial speed is now so high that the axial phase is completed in less than 5 μs and the current is still rising when it is forced down by the radial phase dynamics.

Note 2

A very important point to note in neutron scaling is that there exists some confusion and even misleading information in published literature because of sloppy practice with regards to I_peak and I_pinch. These quantities are sometimes treated as one and the same or when a distinction is attempted there is then confusion between the total current at the time of pinch and I_pinch. For example in the case of PF1000, there appears to be some disappointment (in their publications) that (at 35 kV) with the current at more than 2 MA, Y_n is still at best in the mid 10^11; and not at least an order of magnitude higher that one might expect for currents around 2 MA. However if you numerically run PF1000 at 35 kV you will find that I_pinch is only 1 MA; so we are not surprised that the measured yield is at best an order of magnitude down from what you would expect thinking that your current is around 2 MA. (scaling at Y_n~I⁴ , a factor of 2 in current gives a factor of 16 in the yield; at Y_n~I³ , a factor of 8). So it is important that the thinking of yield should be in terms of I_pinch as the relevant scaling parameter. When using this model code, the distinction of I_pinch and I_peak is clear.

Next we look at the detailed tabulations: As P₀ decreases, I_peak decreases, and continues to decrease, because the increasing axial speed increases the circuit loading, throughout the whole range of pressures. However it is noticed that I_pinch increases from high pressures, peaking in a flat manner at 6 Torr and then decreases sharply as pressure is reduced towards 1 Torr. One factor contributing to the increase is the shift of the pinch time from very late in the discharge (when discharge current has dropped greatly) to earlier in the discharge (when current has dropped less). That is the main factor for I_pinch increasing despite a decreasing I_peak. At low pressures (e.g. 1 Torr), the radial phase now occurs so early that it is forcing the current down early in the discharge. That lowers both the I_peak as well as the I_pinch. These points are clear when you look at the comparative chart of current traces at various pressures.

The radial EINP follows the same pattern as I_pinch, and for the same reasons. The radial EINP computes the cumulative work done by the current sheath (piston) in the radial phases.

Looking at the other quantities, we note that the speeds (axial, radial shock and radial piston) and temperature all continue to rise as pressure lowers; similarly S and maximum induced voltage V also increase as pressure is decreased. Pinch length z_max is almost a constant. Minimum pinch radius and pinch duration continue to decrease; the former due to better compression at higher speeds and the latter due to the increased T. The number density progressively drops, due to the decreasing starting numbers, despite the increasing compression.

From the tabulations of the above numerical experiments, it might be useful to consider the beam-target mechanism which we are using to compute the neutron yield. This is summarized in the following note.

Note 3

(Taken from SP3.doc)

Y_b-t= C_n n_i I_pinch²z_p²(ln(b/r_p))s/V_max^1/2

where s is the D-D fusion cross section. In the range we are considering we may take s~V_maxⁿ where n~2-3; say we take n=2.5; then we have

Y_b-t ~ n_i I_pinch²z_p²(ln(b/r_p)) V_max²

The factor z_p²(ln(b/r_p)) is practically constant.

Thus we note that it is the behaviour of n_{i ,}I_pinch and V_max as pressure changes that determines the way Y_n increases to a maximum and then drops as pressure is changed.

An additional experiment is suggested, in which you can see how numerical experiments on Y_n versus operating pressure compare with measured results in the case of PF400J. This is discussed in Module 7.

4.6 General notes on fitting, yield scaling and applications of the Lee model code

On fitting: In the numerical experiments we soon learn that one is not able to get a perfect fit; in the sense that you can defend it as absolutely the perfect fit. The way to treat it is that one has got a working fit; something to work with; which gives comparable results with experiments; rather than perfect agreement. There is no such thing anyway; experiments in Plasma Focus (i.e. on one PF under consistent conditions) give a range of results; especially in yields (factor of 2-5 range is common). So a working fit should still give results within the range of results of the hardware experiment.

Even though a fit may only be a 'working' fit (as opposed to the hypothetical perfect fit) when one runs a series of well planned numerical experiments one can then see a trend e.g. how properties, including yields, change with pressure or how yields scale with I_pinch, or with L₀ etc. And if carefully carried out, the numerical experiments can provide, much more easily, results just like hardware experiments; with the advantage that after proper reference to existing experiments, then very quickly one can extend to future experiments and predict probable results.

On scaling: Data used for scaling should be taken from yield-optimized (or at least from near optimized) situations. If one takes from the worst case situations e.g. way out in the high pressure or low pressure regions, the yield would be zero for a non-zero I_pinch. Such data would completely distort the scaling picture.

Not only should the pressure be changed, but there should be consideration for e.g. suitable (or even optimized) I_pinch/a; as the value of I_pinch/a would affect the pressure at which optimized S is achieved.

On directions of work and applications: Efforts on the model code may be applied in at least two directions. The first direction is in the further development of the code; e.g. trying to improve the way the code models the reflected shock region or the pinch region.

The second direction is to apply the model to provide a solution to a particular problem. An example was when it was applied to look at expected improvements to the neutron yield of the PF1000 when L₀ is reduced.

Using the model code it was a relatively easy procedure, firing shots as L₀ was reduced in steps; optimizing the various parameters and then looking for the optimized neutron yield at the new value of L₀. When this exercise was carried out in late 2007, for PF1000 at 35 kV, unexpectedly it was found that as L₀ was reduced from 100 nH in steps, in the region around 35 nH, I_pinch achieved a limiting value; in the sense that as L₀ was reduced further towards 5 nH, whilst I_peak continued to increase to above 4 MA, I_pinch dropped slightly from its maximum value of 1.05 MA to just below 1 MA. This Pinch Current Limitation Effect could have considerable impact on the future development of the plasma focus.

On numerical experiments to enhance experience and intuition: Moreover the relationship between I_peak and I_pinch is implicit in the coupling of the equations of circuit and motion within the code which is then able to handle all the subtle interplay of static and dynamic inductances and dynamic resistances and the rapid changes in distributions of various forms of energies within the system. Whilst the intuitive feel of the experienced focus exponents are stretched to the limit trying to figure out isolated or integrated features of these interplays, the simplicity of the underlying physics is captured by the code which then produces in each shot what the results should be; and over a series of shots then reveal the correct trends; provided of course the series is well planned.

So the code may also be useful to provide the numerical experimenter time-compressed experience in plasma focus behaviour; enhanced experience at much reduced time. At the same time the numerical experimenter can in a day fire a number of different machines, without restrictions by time, geography or expense. The problem then becomes one of too much data; sometimes overwhelming the experience and intuition of the numerical experimenter.

On versatility: Your numerical experiments have included examining plasma focus behaviour comparing BIG, medium size and small plasma focus, looking for common and scalable parameters. You studied neutron yields as functions of pressure, comparing computed with experimental data. In 4 modules involving some 12 hours of hands-on work you have ranged over a good sampling of plasma focus machines and plasma focus behaviour.

This was all done with one code the RADPFV5.15de.xls the universal plasma focus laboratory facility. We should have the confidence that if we explore the open experiments suggested in the last module of the advanced course below, that could lead us to new areas and new ideas.

End of Module 4- End of Basic Course in Plasma Focus Numerical Experiments

[Comments and interaction on the course work and other matters related to plasma focus are welcome at anytime]

Reference to this course and the Lee model code should be given as follows:

Lee S. Radiative Dense Plasma Focus Computation Package (2011): RADPF www.plasmafocus.net http://www.intimal.edu.my/school/fas/UFLF/

Also see list of papers at end of Module 1 above

S Lee and S H Saw