and Downloads - Lee model code
· Numerical Experimental Facility
· Simulate any Mathers-type plasma focus, computes dynamics
· Design new plasma focus machines
· Thermodynamics included; H2, D2, Ne, Ar, Xe, He, N2, Kr and D-T
· Model parameters to fit experimental axial, radial phase times
· Radiative phase computes line radiation, recombination and total yield. Computes neutron yield for deuterium operation; based on an improved beam-target model and calibrated at an experimental point of 0.5MA Yn=5.6x10^9 (extended to D-T). Termination of computation of slow phase: transit time for slow disturbance speed. Plasma Self-absorption based on revised equations presented in File 3; appendix by N A D Khattak.
Time guard feature
Choice of Tapered electrode [note: if this feature is not wanted ensure taper? y=1, n=0 cell is filled with a ‘0’]
Quick choice of specified machines; one click loading of chosen machine; at present 3 machines may be click-loaded: the UNU/ICTP PFF, the NX2 and the PF1000 [note: if this feature is not wanted ensure the relevant cell in top RH corner of sheet is filled with a ‘0’]
Acknowledgement: Scroll down to see acknowledgement.
There are altogether 4 files in this package.
File1: PDF File "Description of Radiative Dense Plasma Focus Computation Package": This file
File2: PDF File "Theory of Radiative Plasma Focus Model"
File3: PDF file "Appendix by N A D Khattak".
File7: EXCEL file containing the ACTIVE SHEET AND THE EXECUTABLE (MODIFIABLE) MACRO PROGRAMME CODE.
"Radiative Dense Plasma Focus Computation Code" RADPFV5.15de.c1
In addition, there are files for the computation of thermodynamic data needed for this code.
Hint for downloading the EXCEL FILE: Instead of left click to open the file; it is better to right click and select "save target as"; then choose a suitable location e.g. desktop. The saved EXCEL file will be only about 1M. (see last page for more hints on saving/copying )
These files may also be downloaded from the following URL:
http://www.kirkbyites.net/DPF (only a sketch) or
http://eprints.ictp.it/85/ (containing an earlier version RADPFV5.008)
Files 4 & 5 & 6 contain earlier versions of the code.
A simple 2 phase (axial and radial) model was developed by S.Lee in 1985 as a component of a 3kJ plasma focus experimental package which became known as the UNU/ICTP PFF. This network of basically identical 3kJ PF machines, with different experimental and application emphases, is now operated by groups in countries including Singapore, Malaysia, Thailand, Indonesia, India, Pakistan,, Egypt and Zimbabwe.
The model was written as a 3 phase
(non-radiative) model (in GWBASIC) for an experimental program at the 1991
The present 5-phase package (axial, radial inward shock, radial reflected shock, slow compression radiative and expanded large column phase) is re-written in Microsoft EXCEL VISUAL BASIC in order to make it available for wider usage.
The model may be configured to any conventional Mather-type plasma focus by inputting machine parameters: inductance, capacitance, electrode radii and length. And operating parameters: charging voltage and fill gas pressure. The thermodynamics (specific heat ratio and charge number as functions of temperature) are included for 6 gases namely hydrogen, deuterium, neon, argon and helium. The gases may be selected by simply inputting atomic number, molecular weight and dissociation number (2 for deuterium and hydrogen, 1 for the others).
The results are the following: waveforms for the total discharge current and tube voltage, axial phase trajectory and speed, radial trajectories for the shock front, current sheath and column length and the corresponding speeds, plasma temperature and radiation yields (Bremsstrahlung, line and recombination) and power; and thermodynamic quantities such as specific heat ratios and charge numbers. These are output in graphical as well as tabular forms. Also computed are plasma pinch current and neutron yield, and any distributions, if required.
Are the results any good?
But are there any indications that our computed results are anywhere near the actual results that may be measured on the device in actual operation?
NOT if we just guess the model parameters fm, fc, fmr, fcr; as suggested above. Then the results are just hypothetical; although with experience we may assign some reasonable values of the model parameters for the particular machine in its particular operating conditions. And the results may be valid for planning or designing purposes.
How do we make the results realistic?
The standard practice is to fit the computed total current waveform to an experimentally measured total current waveform.
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.
The exact time profile of the total current trace is governed by the bank parameters namely capacitance Co, external, or static inductance Lo and circuit resistance ro, by the focus tube geometry namely electrode radii, outer ‘b’ and inner anode ‘a’, and the anode length ‘zo’; and on the operational parameters which are the charging voltage Vo and the fill pressure Po and the fill gas. 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. The detailed profile of the discharge current during the pinch phase will also reflect the joule heating and radiative yields. At the end of the pinch phase the total current profile will also reflect the sudden transition of the current flow from a constricted pinch to a large column flow. Thus the discharge current powers all the 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 occurs in the various phases of the plasma focus.
The discharge current of a plasma focus is one of the simplest measurement that is routinely carried out. Yet this easily measured waveform carries all the information of the physical processes in the plasma focus.
Our standard practice for any existing plasma focus is to obtain a measured current trace. Then we fit the computed current trace to the measured current trace. The fitting process involves adjusting the model parameters fm, fc, fmr, fcr one by one, or in combination until the computed current waveform fits the measured current waveform.
Once this fitting is done our experience is that the other computed properties including dynamics, energy distributions and radiation are all realistic
Fitting computed current trace to experimental current trace of existing machine:
The main model parameters are the tube current flow factor CURRF (known to be 0.7 for most machines) and the mass swept-up factor (MASSF, for axial & MASSFR, for radial). These have been pre-selected in the model, but may be adjusted so that the time of focus, and the radial inward shock transit time, fit the experimentally observed times for each machine. The computed current trace may be compared with the experimental current trace.
Features for comparison include current risetime and rising shape, peak current, current 'roll off' and dip, both shape and amplitude. Absolute values should be compared. Our experience with a number of machines shows that the fit is usually very good, occasionally almost exact..
The machine parameters and operating conditions should already have been determined and inputted into the active sheet. The model parameters are then adjusted, one by one, or in combination until best fit is obtained between the computed current trace and the experimental current trace.
Besides the model parameters, sometimes (when all else fails in the fitting process) the inductance (as published or given by the experimenters) needs to be adjusted. The reason is that the inductance Lo may be given as the short circuit bank inductance whereas it should be the ‘static’ inductance of the plasma focus; ie the inductance of the PF before the current sheet moves. Usually also the value of stray resistance ro needs to be guessed at as few experimenters determine this carefully. We usually start with the value of ro as 0.1 of (Lo/Co)0.5; and make small adjustment as necessary. Sometimes, especially for PF’s using very low values of Co, it may also be necessary (when all else fails) to adjust the value of Co (for sub-uF capacitor banks, the closely spaced connecting parallel plates and parallel connecting cables may actually significantly change the value of Co. In one or two cases where there is very good fit in current profiles but the absolute values of currents don’t match, it has been reasonable to suspect that the calibration constant for the current profile has been given wrongly by the experimenter.
Designing a new plasma focus
If a machine has not been built the model may be used to aid design. First use the following rule of thumb procedure [use SI units].
What capacitance ( C ) are you planning?
How low is the inductance (L) you expect to attain?
What maximum voltage (V) do you expect to operate?
Enter these values into the appropriate spaces for these machine parameters.
For the stray (circuit) resistance, take 1/4 the value of (L/C)1/2 .
Estimate the undamped peak current using the formula I=V/(L/C)1/2.
Use (I/a)=250kA max undamped current per cm to assign the value of centre electrode radius 'a'.
Put in double this value for outer electrode radius 'b'.
The length of the electrode may be assigned as 5 times the value of 1.6(LC)1/2 . This length is in cm when the value of (LC)1/2 is expressed in microsecond. (This gives a length which will provide an average axial speed of 5 cm/us which usually gives a peak speed at end of axial phase of 8cm/us. For operation in H2 or D2 20% longer may be better; for neon operation to get suitable line radiation (12-13.5 A) for SXR microlithography purposes, 20% shorter may be better as we require a lower speed to get to the correct level ionization stages. The focus is normally operated so that the start of current dip (signifying the end of the axial phase) occurs at or just after peak current. For argon to generate SXR a high speed is required so use the same length as for H2. For xenon, if the aim is for EUV (around 13 nm) for experiments for NGL the model has predicted a requirement for very low speeds, around 1.3cm/us! So it appears one needs very short anodes, at least 5 times shorter than that needed for D2. However there is not much experimental experience accumulated so far.
For pressure values assign as follows: D2: 4 torr; neon: 1.5 torr; Ar: 0.7 torr.
For xenon, runs on the code suggest several torr to go with the very short anode length.
From the above rule of thumb design parameters, is your PF fat or thin? (according to the ratio length of the centre electrode divided by diameter; for NX2 this ratio is 1.2, fat; for UNU/ICTP PFF this ratio is 17, thin )
If it is fat use the model parameters suggested for the NX2 (These suggested values are tabulated at the top right of the active sheet which appears when you open RADPFV5.008).
If it is thin assign the model parameters closer to the UNU/ICTP PFF.
Run the computation and from results make adjustment to 'a', 'b', length (V may also easily be varied, especially reduced since we have started with max V; C also, use more or less capacitors; careful with L, normally make L as small as possible, but be realistic). Adjust parameters for best results over a range of pressures and gases. Best results could mean strong current dip or biggest line emission in the case of neon, which is useful for developing microlithography SXR sources.
Model may be adapted to suit requirements
The axial phase (trajectory) going into the radial phase (trajectories of shock front, current sheet and length of the focussing column) is portrayed reasonably well. As the radial inward shock goes on-axis, a reflected shock phase follows. The reflected shock moves outward until it hits the incoming radial piston which was moving behind the radial inward shock. Now follows the slow compression, radiative phase.
The radiative phase is the interesting phase, which is presented in the package in a form which gives reasonable results. The slow piston motion is coupled to the rate of change of current, the elongation and power gain/loss due to joule heating, Bremsstrahlung and line emission and recombination losses. Thus radiation collapse (critical current of 1.6 MA for deuterium, but much reduced to possibly below 100kA for neon and argon under certain conditions) is included into the modelling. Reasonable line radiation yields (all lines) are computed for the UNU/ICTP PFF. One may wish to include emission in specific lines. Plasma self-absorption effects (adapted from equations discussed by Khattak in File 3) have been implemented in the code.
There is room for further interesting development. For example, allowing the radiation collapse to couple to a ‘piston’ motion will lead to a huge voltage spike in a ‘high pressure’ regime which is not observed experimentally. In this model, this effect is ‘artificially’ restricted by ‘house keeping’ procedures in the package. It should be looked into further. Instabilities could be introduced into the package by insertion of a suitable ‘anomalous’ resistance time function. This should be coupled into the voltage and current equation; but not into the ‘piston’ equation. The current would rapidly diffuse as this ‘anomalous’ resistance kicks in, causing the necessary abandonment of the concept of a ‘piston’.
The radiative phase is followed by an expanded large column phase, in which the current flows in a large column with the radius of the centre electrode.
The theoretical basis (with all equations) is given in the next (separate) PDF file; File 2., with an appendix by N A D Khattak in PDF file 3.
Most recently a reasonable beam-target mechanism has been incorporated into the code.
This model and code is the result of almost 40 years of experimental work on shock waves and the Plasma Focus. It incorporates the best experience (as understood by me) shared with each and every one of the students and Fellows I have trained and each and every one of the colleagues I have worked with. Some of them are named here: Dr Chow Sai Pew, (the late) Dr Chen You Hor, Yong Yeow Chin, Prof (Dr) Tou Teck Yong, Dr Liu Mahe, Dr Bing Shan, Assoc Prof (Dr) Paul C K Lee, Assoc Prof (Dr) Rajdeep Singh Rawat, Dr Moo Siew Pheng, Assoc Prof (Dr) Kwek Kuan Hiang, Prof (Dr) Wong Chiow San, Dr Chew Ah Chuan, Dr Adrian Serban, Dr Susetyo Mulyodrono, Assoc Prof Zhang Guixin, Prof (Dr) Jalil Ali, Assoc Prof (Dr) Saw Sor Heoh, Dr Alin Patran, Tan Thian Hoo, Dr Suresh Kumar, Chin You Hon, Dr. C.K.Chakrabarty, Ashutosh Srivastava, Assoc Prof (Dr) N A D Khattak.
The package consists of an ACTIVE SHEET in EXCEL linked to a MACRO (where the basic programme is written). The package may be operated from the sheet, without going into the MACRO. The machine parameters may be inputted directly onto the sheet, as may the operating conditions and gas. The model parameters (CURRF, MASSF and MASSFR), if required to be adjusted can also be directly inputted onto the sheet.
After downloading the programme, the EXCEL Sheet appears. The first section, first 19 rows, contain essential information and inputted quantities. Parameters which may be changed directly on the sheet are in bold underlined.
The programme (as downloaded) contains all machine parameters for the UNU/ICTP PFF. Thermodynamic data (for the 3 gases D, Ne and Ar; and xenon) are also preloaded. The pre-selected operational conditions and gas are shown.
The programme may be operated directly from the sheet. Place the cursor in any empty cell; then press “Ctrl + A”
The computation should take less than 1 minute (depending on speed of machine and input parameters).
Results are outputted on the sheet in columns as follows (starting row 20):
A time in microsecond
B current in kA
C tube voltage in kV
D axial position in cm
E axial speed in cm/ms
F time in microsecond (starting on row 20, radial phase)
G time in nanosecond, referenced to start of radial phase
H current in kA (radial phase data only)
I tube voltage in kV (radial phase data only)
J radial shock position, referenced to axis, in mm
K radial piston position, referenced to axis, in mm
L axial position of focus column, referenced to anode end, in mm
M radial shock speed in cm/ms
N radial piston speed in cm/ms
O elongation speed of column in cm/ms
P reflected shock radial position in mm
Q temperature in oK
R Joule heating power in watts
S Bremsstrahlung emission power in watts
T recombination emission power in watts
U line emission power in watts
V sum of S,T & U
W sum of power ie of R & V
X time-integrated Joule heating in Joules
Y time-integrated Bremsstrahlung in Joules
Z time-integrated recombination emission in Joules
AA time-integrated line emission in Joules
AB Total radiation in Joules
AC Total energy gain/loss in J, ie X+AB
AD plasma self-absorption corrected coefficient
AE plasma radiation power (if black body) in W
AF specific heat ratio
AG effective plasma charge number
AH Thermal neutron yield
AI Beam-target produced neutron yield
AJ Total neutron yield
AK plasma ion density per m^3
AM Surface radiation power in W
AN plasma self-absorption corrected coefficient
AO radial phase piston work in J
The number of data rows may go up to 7000. The data is also presented near the top of the sheet in graphical forms. The lines in each figure (identified by colour) plot the following data:
Series 1 dark blue
Series 2 pink
Series 3 yellow
Series 4 light blue
The horizontal axis (for Figs 1 & 2) is time in ms. The other Figures display computed data of the radial phases. Radial phase time scale is in ns, referenced to the start of the radial phase’
Fig 1 (top left)
Fig 2 (top right)
Radial shock position
radial piston position
axial focus length
Radial shock speed
radial piston speed
line radiation energy
Joule heat energy
Line radiation power
Inset Sp Heat Ratio charge number
To get into the code (containing the programme lines)
(from top toolbar) click on 'Tools'
(when menu appears) select 'Macro' and click on it
(when menu appears) select 'Macros' and click on it
(when Macro menu appears) select 'radpf005' by clicking on it
On the right hand panel click on 'step into'
That gets you into the code. You may then modify the code as required.
When finished with the code
click 'red cross' on top right hand corner
click 'OK' when message "This command will stop the Debugger" appears
This gets you back to the active sheet.
The complete computer package (active sheet and macro code) is provided as an EXCEL File RADPFV5.13.2 (latest)
Open the file.
When the EXCEL Sheet appears, press ‘Ctrl + A”. The programme will run and present data in the 8 Figures and inset.
It is recommended that you keep a reference copy of RadPFV5.13.2 (or the current version), which you can refer to. Modified programmes should be kept under another title such as RadPFTrialVersion. (so that you can get back to the original latest version should your trial code not work)
The model has been tested for the UNU/ICTP PFF 3kJ plasma focus, in the following ranges, operating at 14kV.
Deuterium 0.5 to 19 torr
Neon 0.1 to 5.5 torr
Argon 0.1 to 2.5 torr
Any plasma focus will only be able to operate properly within a range of parameters. For example if the parameters are such that the current sheet moves too slowly in the axial phase, by the time the radial phase starts the drive current may have dropped to too low a value and the radial phase cannot complete.
In other words there needs to be a matching between the sum of the characteristic axial & radial times and the characteristic capacitor discharge time.
The code has recently incorporated a Time Match Safeguard. Before axial phase computation is started the code checks that there is suitable matching within a suitable range. This is done by checking that the ratio ALT in the code is LESS than a certain value. We have set the lower limit of ALT at 0.68 for D2 and 0.65 for neon and argon. If too large a pressure is set for selected capacitor voltage and the value of ALT falls below this set value an error message will appear recommending to set the pressure lower or the voltage higher.
If this happens click the 'red cross' on the upper right hand corner, click OK on the 'debugger' message, getting you back to the active sheet and adjust the pressure and or voltage accordingly.
(The value of ALT for each computation is shown in the active sheet.)
Furthermore if you have changed drastically the values of the parameters such as the capacitance C e.g. changing the value of C from 30 uF to 3 uF, you would have reduced the capacitor discharge time by 3 times (square root of 10), hence it may be difficult to match the axial transit time (by simply increasing voltage or reducing pressure) unless you also reduce the anode length by a similar factor. Likewise if you drastically alter the value of the circuit inductance L or the anode radius a you would have to adjust other parameters accordingly. Time matching is crucial for proper computation just as it is in actual operation in the laboratory.
If machines parameters of another machine are entered, keep in the middle (or lower end) of the pressure range. If parameters are inappropriately chosen the Time Match Safeguard will stop the programme execution and give appropriate instructions to remedy.
Another safeguard for inappropriately low operating pressures is also incorporated into the code which stops the execution and gives you appropriate warning.
You may also need to adjust the source data of the figures. Source data has been set for a maximum of 7000 points for Figs 1 & 2, and 6000 points for the other Figures.
The Active sheet comes pre-loaded with parameters of the UNU/ICTP PFF
Operating parameters set as 13kV in neon at 3 torr
If you have difficulties you cannot solve, please e-mail me parameters of your machine.
See the following for a further hint on efficient saving of RADPFV5.08 (applies also to later versions.)
Further hints for efficient saving of RADPFV5.08
3 methods of saving copies: Copy then Paste - small storage required for a copy
Open RADPFV5.008 then Save As- large storage required for a copy
Open RADPFV5.008 then click top Red Cross to exit, click Yes to message 'Save changes…'-small storage required for a copy. If you know why, please let me know; would really appreciate the reasons behind this.
1. The original copy as downloaded (using 'save target as') is around 1M
2 Without opening the file, Use Copy, then Paste to make a copy in a folder separate from your working area. Keep this copy as a reference copy, which you can always return to to make another copy, using Copy then Paste
3. Copy then Paste (without opening the file) will make a copy of 1M.
4. Open the file, then click File, click Save As, will save a copy of several M; with the same content!!
5. If you have opened the file, and made any changes, do not use Save As.
Instead click on red cross at top left hand corner as though to exit. Message 'want to save changes' appears click yes and the changed file will replace the old file keeping storage space to a minimum.
e.g. if you open the file, add a comma to one of the 'unused' cells; click File, click Save As, you will end up with a copy (with an extra inconsequential comma) of several M.
Instead of clicking File, if you click top right hand corner Red Cross, message 'Do you want to save changes…" appears, click yes, the file with the extra comma is saved in place of the old file, without the extra comma and the storage space is still 1M.
Another example: If you open the file, run computation by using CTRL+a, the active sheet now has 8 filled graphs and one filled inset (unlike the original active sheet with all empty graphs). If you save with Save As method you will save a file of perhaps 16M.
If you click on the Red Cross and then Yes to the message you save a file (with same content) of perhaps couple of M. Of course you would no longer have the original file, having replaced it with the file containing the computed results. One more reason to always keep a reserve copy.