Code OK1—Simulation of multi-beam irradiation on a spherical target in heavy ion fusion

Computer Physics Communications 157 (2004) 160–172 www.elsevier.com/locate/cpc

Code OK1—Simulation of multi-beam irradiation on a spherical target in heavy ion fusion ✩ A.I. Ogoyski a,b,∗ , T. Someya a , S. Kawata a a Department of Energy and Environmental Science, Utsunomiya University, Utsunomiya 321-8585, Japan b Department of Physics, Varna Technical University, Varna 9010, Bulgaria

Received 7 June 2003; received in revised form 26 August 2003; accepted 19 October 2003

Abstract Code OK1 is a fast and precise three-dimensional computer program designed for simulations of heavy ion beam (HIB) irradiation on a direct-driven spherical fuel pellet in heavy ion fusion (HIF). OK1 provides computational capabilities of a three-dimensional energy deposition profile on a spherical fuel pellet and the HIB irradiation non-uniformity evaluation, which are valuables for optimizations of the beam parameters and the fuel pellet structure, as well for further HIF experiment design. The code is open and complete, and can be easily modified or adapted for users’ purposes in this field. Program summary Title of program: OK1 Catalogue identifier: ADST Program summary URL: http://cpc.cs.qub.ac.uk/summaries/ADST Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland Computer: PC (Pentium 4, ∼1 GHz or more recommended) Operating system: Windows or UNIX Program language used: C++ Memory required to execute with typical data: 911 MB No. of bits in a word: 32 No. of processors used: 1 CPU Has the code been vectorized or parallelized: No No. of bytes in distributed program, including test data: 16 557 Distribution format: tar gzip file Keywords: Heavy ion beam, inertial confinement fusion, energy deposition, fuel pellet Nature of physical problem: Nuclear fusion energy may have attractive features as one of our human energy resources. In this paper we focus on heavy ion inertial confinement fusion (HIF). Due to a favorable energy deposition behavior of heavy ions in matter [J.J. Barnard et al., UCRL-LR-108095, 1991; C. Deutsch et al., J. Plasma Fusion Res. 77 (2001) 33; T. Someya et al.,

✩ This paper and its associated computer program are available via the Computer Physics Communications homepage on ScienceDirect (http://www.sciencedirect.com/science/journal/00104655). * Corresponding author. E-mail address: [email protected] (S. Kawata).

0010-4655/$ – see front matter  2003 Elsevier B.V. All rights reserved. doi:10.1016/S0010-4655(03)00492-2

A.I. Ogoyski et al. / Computer Physics Communications 157 (2004) 160–172

161

Fusion Sci. Tech. (2003), submitted] it is expected that heavy ion beam (HIB) would be one of energy driver candidates to operate a future inertial confinement fusion power plant. For a successful fuel ignition and fusion energy release, a stringent requirement is imposed on the HIB irradiation non-uniformity, which should be less than a few percent [T. Someya et al., Fusion Sci. Tech. (2003), submitted; M.H. Emery et al., Phys. Rev. Lett. 48 (1982) 253; S. Kawata et al., J. Phys. Soc. Jpn. 53 (1984) 3416]. In order to meet this requirement we need to evaluate the non-uniformity of a realistic HIB irradiation and energy deposition pattern. The HIB irradiation and non-uniformity evaluations are sophisticated and difficult to calculate analytically. Based on our code one can numerically obtain a three-dimensional profile of energy deposition and evaluate the HIB irradiation non-uniformity onto a spherical target for a specific HIB parameter value set in HIF. Method of solution: OK1 code is based on the stopping power of ions in matter [J.J. Barnard et al., UCRL-LR-108095, 1991; C. Deutsch et al., J. Plasma Fusion Res. 77 (2001) 33; T. Someya et al., Fusion Sci. Tech. (2003), submitted; M.H. Emery et al., Phys. Rev. Lett. 48 (1982) 253; S. Kawata et al., J. Phys. Soc. Jpn. 53 (1984) 3416; T. Mehlhorn, SAND80-0038, 1980; H.H. Andersen, J.F. Ziegler, Pergamon Press, 1977, p. 3]. The code simulates a multi-beam irradiation, obtains the 3D energy deposition profile of the fuel pellet and evaluates the deposition non-uniformity. Restrictions on the complexity of the problem: No Typical running time: The execution time depends on the number of beams in the simulated irradiation and its characteristics (beam radius on the pellet surface, beam subdivision, projectile particle energy and so on). In almost of the practical running tests performed, the typical running time for one beam deposition is less than 2 s on a PC with a CPU of Pentium 4, 2.2 GHz (e.g., in Test 2 when the number of beams is 600, the running time is about 18 minutes). Unusual features of the program: No  2003 Elsevier B.V. All rights reserved. PACS: 52.58.Hm; 28.52.Cx; 29.27.Eg Keywords: Heavy ion beam; Inertial confinement fusion; Energy deposition; Fuel pellet

1. Introduction Nuclear fusion energy may have attractive features as one of our human energy resources. In the nuclear fusion the fuel is consisted of isotopes of hydrogen, for example, deuterium or tritium. The fuel exists globally in the sea or can be bred. Once the fusion energy to become available for humanity, the energy problem we will meet in the near future will be solved. At present schemes of magnetically confinement fusion and inertial confinement fusion (ICF) have been studied intensively. In this paper we focus on heavy ion ICF. Due to a favorable energy deposition behavior of heavy ions in matter, high accelerator efficiency and recent research explorations [1–3], it is expected that HIB would be one of energy driver candidates to operate a future inertial confinement fusion power plant. Preferable features in HIF include a direct couple of the HIB energy to a fuel pellet and a relatively deep-area heating of an energy deposition layer. For a successful fuel ignition and fusion energy release, a stringent requirement is imposed on the HIB irradiation non-uniformity, which should be less than a few percent [3–5]. In order to meet that requirement we need to evaluate the non-uniformity of realistic HIB irradiation and energy deposition pattern onto a target. In HIF, HIBs would be expected to be neutralized in a HIF reactor chamber and HIB ion trajectories may be calculated ballistically. Then Each HIB is focused onto a small target, located at the reactor center. In this paper we present our code for HIB irradiation computations in HIF. The computer code is oriented for HIB irradiation energy deposition and non-uniformity evaluation for a direct-driven spherical pellet. HIB ions impinge the pellet surface and deposit their energy in a relatively deep and wide area. Therefore the HIB irradiation non-uniformity should be evaluated in the volume of the deposition area in the energy absorber layer. The stopping power is solved for HIB ions by including target bound and free electrons, target nuclei and target ions. In our code each HIB is divided into many beamlets and each beamlet energy deposition is computed precisely. The total HIB irradiation pattern is summarized over all the HIBs irradiated onto a target and the HIB irradiation nonuniformity is evaluated by the total relative root-mean-square (RMS) deviation, and total relative peak-to-valley

162

A.I. Ogoyski et al. / Computer Physics Communications 157 (2004) 160–172

(PTV) deviation. Based on our code one can obtain a three-dimensional picture of energy deposition and evaluate the HIB irradiation non-uniformity onto a spherical target for a specific HIB parameter value set in HIF. The existing irradiation systems can be optimized and farther HIF experiments can be designed.

2. OK1 algorithm description The code OK1 is based on a fast and precise algorithm. Its principal steps are: (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15) (16) (17)

Initialize the beam and target parameters. Take the first bunch. Set the target temperature. Take the first beam from the fixed irradiation system. Divide the beam into 316 beamlets (see Fig. 4). Take the external pellet shell. Find the beamlet center crossing points with this pellet shell (see Fig. 6). Calculate the energy density deposition and fill the nearest pellet mesh with it. Fill the meshes in the beamlet area with the same energy density (see Fig. 7). Take the next pellet shell (till the completely ion energy deposition) and go to step 7. Take the next beamlet (till the last one) and go to step 6. Correct the one-beam energy deposition systematic error. Take the next beam (till the last one) and go to step 5. Check up the conservation of energy. Take the next bunch (till the last one) and go to step 3. Save the results. End.

3. Included files 1. StoppingPower1.cpp: The file contains a function Stop1. This function serves a heart of the OK1 code and describes the energy deposition model. It calculates the stopping power from the projectile ions into the solid target. The one-ion stopping power is considered to be a sum of the deposition energy in target nuclei, target bound and free electrons and target ions [6]. EStPower = Enuc + Ebound + Efree + Eion .

(1)

The nuclear stopping power process is effective in the end of the stopping range and describes the elastic Coulomb collisions between the projectile ions and target nuclei [7,8]: √ 
Enuc = 10−7 Cn1 E exp −45.2(Cn2E)0.277 , (2) 3/2  1/2 
2/3 Zb Zt Ab 2/3 −3/2 Cn1 = 4.14 × 106 ρ (3) Zb + Zt , Ab + At Ab  
2/3 Ab At 3/2 2/3 −1/2 −6 (Zb Zt )−1 Zb + Zt , Cn2 = 10 ρ (4) Ab + At where Ab and At are the projectile and target atomic weights, Zb and Zt —the projectile and target atomic numbers, E—the projectile ion energy and ρ—the target density.

A.I. Ogoyski et al. / Computer Physics Communications 157 (2004) 160–172

163

The Linhard and Bethe–Bloch equations describe the bound electron stopping power process. The model uses the Thomas–Fermi description of the electron clouds and includes the Coulomb collisions as well as excitation and ionization processes [7–10]: Ebound = ELSS ∪ EBethe ,

1/3

with LSS validity criterion: Zb


137

Vbeam , c

(5)

7/6

Zb 2e2 aB neb Vbeam (aB and VB are the Bohr radius and velocity), 2/3 ε0 VB (Z + Zt2/3 )3/2 b   2 V2 2me Vbeam e4 neb Zeff 2b − ln(1 − β) − β − δ , β = beam ln . EBethe = 2 2 Ui

c2 4πε0 me Vbeam

ELSS =

The effective projectile charge is    1 Vbeam Zeff b = Zb 1 − 1.034 exp −137.04 0.69 c Zb

(6) (7)

(8)

by the Brown and Moak expression [6–12]. The shell corrections coefficients δ and the average ionization potential Ui of the target are taken from Andersen and Ziegler [7]. The temperature dependence of target ionization Zt (T ) is obtained solving the equation of state (EOS) based on the Thomas–Fermi model [13]. The free electron stopping power plays an important role with the rise of the target temperature [14,15]. e4 nef Zeff 2b (Ge Le + He ), 4πε02 me V re2      2 Xe X2 the Chandrasekhar function, Xe exp − e Ge = erf √ − π 2 2     λD Xe3 Xe2 X4 exp − Le = ln , He = − √ + 4 e . bmin 2 Xe + 12 3 2π ln Xe Efree =

(9) Xe =

V re , Ve

(10) (11)

Here nef is the free electron number density, V re is the relativistic projectile velocity, Ve is the electron thermal velocity, bmin is the distance of closest approach and λD is the Debye length. The free ion stopping power is described analogically by Eqs. (12)–(14). e4 na Zeff 2b Zt2 (Gi Li + Hi ), 2 4πε02 Mt Vbeam      Xi2 2 Xi Vbeam , Xi = Xi exp − Gi = erf √ − , π 2 Vi 2     X4 X3 X2 λD , Hi = − √ i Li = ln exp − i + 4 i , bmin 2 Xi + 12 3 2π ln Xe Eion =

(12) (13) (14)

Mt is the target atomic mass, na is the target atom number density and Vi is the ion thermal velocity. The stopping power for one Pb projectile ion in an Al target in a temperature range of 1 eV–1 KeV is plotted in Fig. 1. 2. InputOK1.h: The input data file contains: projectile and target matter parameters, HIB parameters, pellet mesh parameters, pellet structure and reactor chamber parameters. 3. InitOK1.h: The file initializes the OK1 code and function Stop1. It contains the necessary physical constants, procedure declarations and so on. 4. HIFScheme.h: The file contains all direction coordinates of the irradiation schemes used in our simulations. The beam directions are defined in a spherical coordinate system linked to the reactor chamber center (CS in Fig. 5). 5. Derf.c: The file contains the error function in double precision.

164

A.I. Ogoyski et al. / Computer Physics Communications 157 (2004) 160–172

Fig. 1. Stopping power for one Pb projectile ion in an Al target.

Fig. 2. One-beamlet trajectory.

4. OK1 procedures description 1. Procedure Irradiation(): This procedure organizes a one-beam energy deposition process. 2. Procedure InitEdp1(): This procedure calculates one-beamlet kinetic energy for each bunch. Eprt is the energy of one projectile ion, nE is a particle number parameter, nBunch is the particle number in one bunch and tdM is a beam emittance parameter. 3. Procedure InitE(): This procedure initializes a one-beam energy deposition array Eij k . 4. Procedure InitET(): This procedure initializes a total energy deposition array ET ij k . 5. Procedure InitMesh(): This procedure initializes the coordinates of each mesh cell and store its in arrays rMeshi , ThetaMeshj and PhiMeshk . 6. Procedure InitSurface(): This procedure calculates the arc length rS between the beam center and each one of beamlet crossing points with the pellet surface, solving the system of Eqs. (15).   rS = αRplt ,     y = R sin α, plt (15)  x = R plt cos α,     f +Rplt = f +x . Rbm

y

Rplt is the radius of a fixed pellet surface, f is the focal point distance, Rbm is a beamlet distance from the beam center and rS is the seeking arc length (see Fig. 2).

A.I. Ogoyski et al. / Computer Physics Communications 157 (2004) 160–172

165

7. Procedure InitVMesh(): This procedure calculates the volume of each mesh cell and stores it in V Mij k array, solving the system of Eqs. (16).  V in r = rMesh − dr ,  i i  2   V ex r = rMesh + dr ,   i i  2 (16) , θ1j = j dθ 2   dθ    θ2j = (j +
1) 2 ,  
  4π V ex 3 ri − V in ri3 sin2 θ2j − sin2 θ1j . Vij = 3φ m 8. Procedure ETotal(): This procedure adds the energy deposited from each beam in the total energy deposition array. 9. Procedure InitBeam(): This procedure initializes the beam radius and the beam coefficients A1 , A2 , A3 , A4 and A5 , which serve for the Gauss distribution calculation of the beam particles. The beamlet cell traversal area S1 and volume V1 are also calculated. 10. Procedure Focus(): This procedure calculates the focal point distance f as follows: f = (Rb Rch − R0 Ren )/(Ren − Rb ), where Rb is the beam radius on a tangential α-plane, Rch is the reactor chamber radius, R0 is the external pellet radius and Ren is the beam radius on the chamber entrance (see Fig. 3). 11. Procedure fDis(): This procedure calculates the beam particle number density distribution coefficients Fij (see Fig. 4). The beam density distribution has not axial symmetry in general. If fD (x) and fD (y) are the density distribution functions in transversal horizontal and transversal vertical directions, j +1 i+1

fD (x) dx i fD (y) dy j Fij = 1. , with normalizing condition Fij = R (17) R 4 0 b fD (x) dx 0 b fD (y) dy i,j

Fig. 3. One-beam irradiation onto a pellet with a HIB divergence.

Fig. 4. One-beam division to many beamlets.

166

A.I. Ogoyski et al. / Computer Physics Communications 157 (2004) 160–172

In case of Gauss distribution and axial symmetric beam equation (17) takes an explicit form: xj2 Fij =

√C xj1 s 2π


x 2  yi2 C
y2  √ exp − 2s exp − 2s 2 dx yi1 2 dy s 2π ,
Rb C
x 2  2 4 0 √ exp − 2s 2 dx

(18)

s 2π

xj 1 = bb(mb − 2j ), yi1 = bb(mb − 2i),

xj 2 = bb(mb − 2j + 2),

(19)

yi2 = bb(mb − 2i + 2),

where s is a standard deviation, Rb is the beam radius, mb is a beam divider parameter. The normalizing coefficient C makes the denominator of Eq. (18) equal to 1 (e.g., for Rb = 3.3 mm, mb = 20 and s = 0.55Rb , C takes value of 1.074). In the code Eq. (18) is solved as follows:    Fij = A5 Fxj Fyi ,   

 Fxj = A1 xj 2 − xj 1 − A2 xj32 − xj31 − A3 xj52 − xj51 + A4 xj72 − xj71 , (20)    F = A y − y − A y 3 − y 3 − A
y 5 − y 5  + A
y 7 − y 7 , yi 1 i2 i1 2 i2 3 i2 4 i2 i1 i1 i1 a2

a2 where A1 = √C , A2 = a32 , A3 = 3.3333 , A4 = 142 and a2 = 2s12 . The coefficient A5 compensates the beam s 2π discretization error and assures the normalizing condition of Eq. (17). 12. Procedure Divider(): This procedure organizes the beam division on beamlets. Each beamlet deposits its energy in the target independently. 13. Procedure kBunch(): This procedure calculates the energy deposition coefficients Lij and Kij for each fixed mesh cell. Lij is the total trace length from one-particle traversal plane of each fixed beamlet, dij is the effective particle diameter and Kij is the number of traversal planes for that beamlet. 14. Procedure PointC(): We use two spherical coordinate systems, one linked to the reactor chamber center— Chamber System (CS) and another linked to the pellet center—Pellet System (PS). Each beam axis is fixed by the coordinates θbeam and φbeam in CS (see Fig. 5). This procedure calculates the coordinates of the beam center on the pellet surface rC , θC and φC in PS. 15. Procedure PointF(): This procedure calculates the coordinates of the focal point F : rf , θf and φf in PS (see Fig. 5).

Fig. 5. The beamlets coordinate schemes.

A.I. Ogoyski et al. / Computer Physics Communications 157 (2004) 160–172

167

16. Procedure PointAlpha(): We fix each beamlet trajectory by two points in PS—the crossing point with a tangential α-plane and the common focal point F (see Fig. 5). This procedure calculates the coordinates (r, θ, φ) to each α-point in PS. 17. Procedure CrossP(): A beamlet trajectory is fixed by point A(r, θ, φ) and point F (rf , θf , φf ) in PS, and crosses the pellet surface in two points—B and C (see Fig. 6). This procedure finds the coordinates of the near crossing point B(rCP , θCP , φCP ) for each beamlet, by solving the system of Eq. (21). rCP sin θCP cos φCP − r sin θ cos φ rCP sin θCP sin φCP − r sin θ sin φ rCP cos θCP − r cos θ = = . rf sin θf cos φf − r sin θ cos φ rf sin θf sin φf − r sin θ sin φ rf cos θf − r cos θ

(21)

18. Procedure NearP(): This procedure finds the coordinates of the nearest mesh cell center (rNear , θNear , φNear ) for each crossing point B(rCP , θCP , φCP ). 19. Procedure BeamletEd(): This procedure calculates a one-beamlet energy density deposition in the pellet meshes. 20. Procedure rDis(): This procedure calculates the distance between two points. 21. Procedure Edep1(): This procedure fixes the pellet structure. The code is made for a multi-structure pellet. In our tests we employ an Al monolayer structure, defined in the main procedure. 22. Procedure jkSp(): Each beamlet has the same traversal area S1 but impinges on a different mesh grid density (see Fig. 7). This procedure calculates the spread coefficients jSp in θ -direction and kSp in φ-direction for each beamlet across each one of deposition shells. 23. Procedure Spread(): This procedure chooses the neighbor mesh cells which should be filled with the same energy density like the beamlet center cell (see Fig. 7). 24. Procedure ESp(): This procedure fills the chosen mesh cells (by procedure Spred()) with the calculated energy density. 25. Procedure SD(): This procedure calculates the total relative root-mean-square (RMS) deviation and the total relative peak-to-valley (PTV) deviation as follows:   nθ nφ nr 2

1  Ei k ( E i − Eij k )  j wi σRMS i , σRMS i = , wi = , σRMS = (22) E i nθ nφ E i

σPTV =

nr

i

wi σPTV i ,

σPTV i =

Eimax − Eimin , 2 E i

(23)

where σRMS i is the relative RMS deviation on the ith spherical shell layer of the fuel pellet, Ei is the energy deposited on that shell and E is the total energy deposition, σPTV i is the relative PTV deviation on the ith shell layer, Eimax and Eimin are the maximal and minimal energy deposited on that shell. The energy deposition layer

Fig. 6. One-beamlet crossing points in the pellet.

Fig. 7. One-beamlet energy deposition in a pellet shell.

168

A.I. Ogoyski et al. / Computer Physics Communications 157 (2004) 160–172

is divided into a mesh structure in the r, θ and φ directions with the following mesh numbers: nr = 50, nθ = 90, nφ = 90. Each shell layer has a thickness dr of 0.02 mm. 26. Procedure ECL(): This procedure checks up the conservation of energy. The total energy deposition is compared with the total energy of the impinging beams. The error is controlled to be less than 0.01%. 27. Procedure ADJ(): This procedure corrects the code systematic error due to the relative pellet mesh roughness. The correction includes also special cases when the beamlet touches the pellet surface and only a part of its particles impinge on the pellet layer. After the execution of that procedure the systematic error is reduced to less than 0.1%. 28. Procedure Hole(): In an exceptional deposition condition case it may be possible to have holes in which no beam energy deposition appears in one or few meshes inside of one-beam deposition area, because of the beamdivision roughness. This procedure finds the holes and fills them with a suitable interpolated energy. 29. Procedure SPL(): This procedure makes a surface interpolation smoothing for procedure Hole(). 30. Procedure SPLr(): This procedure makes a radial interpolation smoothing for procedure Hole(). 31. Procedures jkmM(), Espl(), CheckE(), CheckR(): These procedures help procedure Hole() to perform only in the beam deposition area. 32. Procedure View3D(): This procedure transfer the calculated energy deposition value in each pellet cell into that in the Cartesian PS coordinates, for the convenience of further visualization processes.

5. Instructions for the user To run the OK1 code under conditions different from those of the tests described here, a user can control the simulation in 5 principal ways: 1. To specify the HIB parameters and HIB irradiation system. (a) Projectile ion type: Three projectile ion types are included in OK1—Pb, U and Cs. The user can choose one of them or add other species expanding the arrays aZb and aAb. (b) Beam radius at the reactor chamber entrance can be specified by parameter Ren . (c) Beam subdivision in beamlets can be specified by the parameter mb whose value should be an integer even number greater than 4. (d) Beam subdivision in bunches: Each HIB is subdivided in 10 bunches in longitudinal with the length of 83.7 mm and particle numbers stored in the array nBunch. This prescription is prepared to simulate the HIB total energy deposition profile and non-uniformity including the HIB pulse shape. These parameters are obtained from the initial HIB characteristics: the maximal current is 5 kA, the pulse width is 10 ns and the rise and fall times are 2 ns. Users can specify other HIB parameters and bunch subdivision ways changing the arrays lBunch and nBunch. (e) Beam particle number density distribution can be selected by parameter bndd (K–V or Gauss distribution). The user can specify other beam number density distribution changing the coefficients A1 , A2 , A3 , A4 and A5 as it is described in Section 4.11. (f) Beam particle energy distribution is the Maxwell one with the mean particle energy Eb of 8 GeV, minimal energy of 0.6 × Eb and maximal energy of 1.5 × Eb . The beam energy interval is subdivided in 10 intervals and particle number coefficients and energies in each interval are stored in the arrays MxwN and MxwE for two HIB temperatures—100 and 300 MeV in the present OK1 as typical example temperatures for Pb HIBs. The user can select other HIB temperature or specify another HIB energy distributions adding the new values of the arrays. (g) Beam emittance: To specify the HIB emittance, the user should choose first the mean beam radius Rb on the α-plane, calculate the divergence angle αdvr by Eq. (24) and find the minimal and maximal beam radii Rb corresponding to the focal point distances f1 and f2 (see Fig. 3). After that, the user should subdivide each

A.I. Ogoyski et al. / Computer Physics Communications 157 (2004) 160–172

169

HIB by focal point distance into tdM beams and impinge the pellet with the all (see Test 2). π εr = (24) αdvr Ren [mm-mrad]. 0.18 (h) HIB irradiation schemes: The file HIFScheme.h contains 12, 20, 32, 60 and 120-beam irradiation schemes. The user can choose one of them or add other HIB irradiation schemes supplementing that file. 2. To specify the reactor chamber radius. The user can specify the chamber radius changing the parameter Rch . 3. To specify the pellet mesh subdivision. The pellet mesh numbers in r, θ, φ directions in PS are rM, ThetaM and PhiM with the steps dr, dTheta and dPhi, respectively. The user can specify other pellet mesh subdivision ways. 4. To specify the target pellet structure: OK1 includes an example Pb–Al-foam Al–Al–DT structure [3] defined by parameters R0, R1, R2, R3, R4 and R5. Changing these parameters the user can specify other pellet structures. A simple Al mono-structure is used for the tests. The user can add other target matters expanding the arrays aZt0, aZtm, aAt, aUi, aro and SC. The shell corrections coefficients (in the array SC) and the average ionization potential (in the array aUi) of the target are taken from Andersen and Ziegler [7]. 5. To specify the pellet displacement from the reactor chamber center. The parameters dx, dy and dz fix the pellet displacement from the chamber center in the Cartesian PS coordinates. All desirable changes should be made in file InputOK1.h (except for the beam number density distribution and irradiation system).

6. Testing the program OK1 6.1. Test 1 This test illustrates a one-beam energy deposition calculation. The beam has following parameters: the beam radius in the entrance of the chamber Ren = 35 mm; the beam radius on α-plane Rb = 3 mm; the beam particle density distribution is in the Gaussian and all projectile ions have the same energy of 8 GeV (a cold beam). The

Fig. 8. A cross-section diagram of one-beam energy deposition profile in an Al layer.

170

A.I. Ogoyski et al. / Computer Physics Communications 157 (2004) 160–172

Table 1 Output Test 1 (part of the file View3D.dat) . . . . . . . . . X . . . .. . . . . . . . . . . . . Y . . . . . . . . . . . . . . . Z . . . . . . . . . .. dE/dV [J/mm3 ] .............................................................................. −3.447633e−001 2.693587e−001 3.563241e+000 6.855623e+001 −3.627130e−001 2.446530e−001 3.563241e+000 7.777765e+001 −3.788956e−001 2.187555e−001 3.563241e+000 8.680845e+001 −3.932322e−001 1.917922e−001 3.563241e+000 9.491071e+001 −4.056531e−001 1.638945e−001 3.563241e+000 1.961330e+002 −4.160976e−001 1.351983e−001 3.563241e+000 2.016776e+002 −4.245150e−001 1.058435e−001 3.563241e+000 2.061795e+002 −4.308641e−001 7.597301e−002 3.563241e+000 2.086139e+002 −4.351142e−001 4.573238e−002 3.563241e+000 2.131519e+002 −4.372444e−001 1.526895e−002 3.563241e+000 2.168062e+002 −4.372444e−001 −1.526887e−002 3.563241e+000 2.208562e+002 −4.351142e−001 −4.573230e−002 3.563241e+000 2.200691e+002 −4.308642e−001 −7.597293e−002 3.563241e+000 2.173075e+002 ..............................................................................

Fig. 9. σRMS and σPTV versus the pellet radius and the deposition energy profile for the 12-beam irradiation system.

beam is divided into 316 beamlets (mb = 20). The energy deposition proceeds in an Al monolayer pellet with a thickness of 1 mm and an external radius of 4 mm. The results are transferred in the Cartesian PS coordinates and saved by the procedure View3D() in the file View3D.dat (23 MB) and visualized in Fig. 8. The typical running time for one beam deposition is less than 2 s. The procedure View3D() transfers and saves the results for about 6 s. To run the Test 1, the area “Test 1” should be enabled and “Test 2” disabled in the MAIN procedure. A part of the output file View3D.dat is listed in Table 1. 6.2. Test 2 This test illustrates a 12-beam energy deposition calculation. Each beam has the mean radius on α-plane Rb = 3.7 mm. The beam emittance at the reactor chamber entrance is εr = 3.5 mm-mrad with a divergence angle αdvr = 5.74 × 10−3 deg. (see Fig. 3). The minimal beam radius is 3.5 mm and the maximal is 3.9 mm. The beam particle density distribution in one beam is in the Gaussian. The projectile ion energy is in the Maxwell distribution with a beam temperature of 300 MeV. Including the Maxwell distribution (see Section 5.1(f)) and the beam emittance (see Section 5.1(g)), each beam is divided additionally into 10 beams in the energy space and into 5 beams by the focal point distance. Thus the number of beams is 600 and the number of beamlets becomes 189 600 (i.e. 12 × 10 × 5 × 316). The running time for Test 2 is about 18 minutes. The results are plotted in Fig. 9

A.I. Ogoyski et al. / Computer Physics Communications 157 (2004) 160–172

171

Table 2 Output Test 2 (file SD12.dat) Edpi [J]

σRMS i [%]

σPTV i [%]

33351.820813 34477.814218 35631.632053 36808.808394 37980.331521 39063.444640 39905.722204 40364.135298 40288.505904 39575.380785 38169.934474 36067.851079 33459.017240 30739.514196 27446.920367 24256.037037 21418.197024 18378.720788 15797.503110 13110.225393 10962.480826 8707.520981 7022.916194 5287.794297 4260.235723 3191.869731 2134.985982 1552.633553 997.532796 492.510839 224.094315

3.866006 4.175388 4.367818 4.493307 4.469179 4.378458 4.175159 3.732082 3.252660 2.847215 2.600424 2.503364 2.425709 2.460260 2.583174 2.704253 2.773237 2.923529 2.965651 3.068203 3.076494 3.079842 3.520733 4.125071 7.790947 11.893013 14.869986 23.760878 33.166539 37.817662 43.214892

8.476398 9.119073 9.487664 9.722140 9.500493 9.377880 8.934259 7.840950 6.755600 5.676428 5.024234 4.863256 4.663019 4.962266 5.243288 5.686210 6.250784 6.793264 7.395218 7.778643 7.900808 8.458377 8.117253 10.456936 17.489297 24.387743 30.383815 46.419335 65.607871 74.681437 84.956350

σRMS 3.65

σPTV 7.79

errors 0

and saved in file SD12.dat listed in Table 2. The first file column contains the energy deposited in each pellet shell Edp[J], the second σRMS i [%] and the third σPTV i [%]. To run the Test 2 the area “Test 2” should be enabled and “Test 1” disabled in the MAIN procedure. The both tests were performed on C++ on Windows (the PC processor is Pentium 4, 2.2 GHz).

Acknowledgements This research is partly supported by the Japan Society for the Promotion of Science. The authors also would like to present our thanks to all the friends in a HIF research group in VNL, USA, a HIF research group in Japan and Prof. A. Blagoev for their fruitful discussions on this subject.

172

A.I. Ogoyski et al. / Computer Physics Communications 157 (2004) 160–172

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

J.J. Barnard, et al., Lawrence Livermore National Laboratory Research Report, UCRL-LR-108095, 1991. C. Deutsch, et al., J. Plasma Fusion Res. 77 (2001) 33. T. Someya, et al., Fusion Sci. Tech. (2003), submitted for publication. M.H. Emery, et al., Phys. Rev. Lett. 48 (1982) 253. S. Kawata, et al., J. Phys. Soc. Jpn. 53 (1984) 3416. T. Mehlhorn, Sandia Report, SAND80-0038, 1980. H.H. Andersen, J.F. Ziegler, Pergamon Press, 1977, p. 3. P. Wang, et al., Phys. Plasmas 5 (1998) 2977. F.C. Young, et al., Phys. Rev. Lett. 49 (1982) 549. K.A. Long, N.A. Tahir, Phys. Lett. A 91 (1982) 451. K.A. Long, N.A. Tahir, Phys. Rev. A 35 (1987) 2631. J.A. Swegle, Sandia Report, SAND82-0072, 1982. A.R. Bell, Rutherford Laboratory Report, RL-80-091, 1981. T. Peter, J. Meyer-ter-Vehn, Phys. Rev. A 43 (1991) 1998. T. Peter, J. Meyer-ter-Vehn, Phys. Rev. A 43 (1991) 2015.