Routine to calculate parameters of a heavy ion driver suitable for inertial fusion energy author: P J Knight, CCFE, Culham Science Centre aaion : input real : Ion mass (amu) bmax : input real : Maximum field at the superconductor (T) dpp : input real : Beam momentum spread dtheta : input real : Pointing error (rad) edrive : input real : Driver energy (J) emitt : input real : Normalised beam emittance (metre-rad) etai : input real : Axial quadrupole packing fraction at injection lf : input real : Distance from final focussing quad to target (m) nbeams : input integer : Number of beams qion : input real : Ion charge state sigma : input real : Depressed tune (incl. space charge effects) (deg) sigma0 : input real : Phase advance per lattice period (tune) (deg) tauf : input real : Post-acceleration pulse duration (s) theta : input real : Final focussing half angle (rad) vi : input real : Injection voltage (V) gain : output real : Target gain etadrv : output real : Driver efficiency This routine calculates the parameters of a heavy ion driver suitable for inertial fusion energy.
Currently, the complicated model taken from the reference is not
complete, so it is recommended that the simple model is used
(set ISIMP=1
)
Heavy-ion Driver Design and Scaling, R. Bieri et al.,
Fusion Technology, vol.21 (1992) 1583
Meier and Bieri, Fusion Technology, vol.21 (1992) 1547
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=dp), | intent(in) | :: | aaion | |||
real(kind=dp), | intent(in) | :: | bmax | |||
real(kind=dp), | intent(in) | :: | dpp | |||
real(kind=dp), | intent(in) | :: | dtheta | |||
real(kind=dp), | intent(in) | :: | edrive | |||
real(kind=dp), | intent(in) | :: | emitt | |||
real(kind=dp), | intent(in) | :: | etai | |||
real(kind=dp), | intent(in) | :: | lf | |||
integer, | intent(in) | :: | nbeams | |||
real(kind=dp), | intent(in) | :: | qion | |||
real(kind=dp), | intent(in) | :: | sigma | |||
real(kind=dp), | intent(in) | :: | sigma0 | |||
real(kind=dp), | intent(in) | :: | tauf | |||
real(kind=dp), | intent(in) | :: | theta | |||
real(kind=dp), | intent(in) | :: | vi | |||
real(kind=dp), | intent(out) | :: | gain | |||
real(kind=dp), | intent(out) | :: | etadrv |
subroutine iondrv(aaion,bmax,dpp,dtheta,edrive,emitt,etai,lf, &
nbeams,qion,sigma,sigma0,tauf,theta,vi,gain,etadrv)
!! Routine to calculate parameters of a heavy ion driver
!! suitable for inertial fusion energy
!! author: P J Knight, CCFE, Culham Science Centre
!! aaion : input real : Ion mass (amu)
!! bmax : input real : Maximum field at the superconductor (T)
!! dpp : input real : Beam momentum spread
!! dtheta : input real : Pointing error (rad)
!! edrive : input real : Driver energy (J)
!! emitt : input real : Normalised beam emittance (metre-rad)
!! etai : input real : Axial quadrupole packing fraction at injection
!! lf : input real : Distance from final focussing quad to target (m)
!! nbeams : input integer : Number of beams
!! qion : input real : Ion charge state
!! sigma : input real : Depressed tune (incl. space charge effects) (deg)
!! sigma0 : input real : Phase advance per lattice period (tune) (deg)
!! tauf : input real : Post-acceleration pulse duration (s)
!! theta : input real : Final focussing half angle (rad)
!! vi : input real : Injection voltage (V)
!! gain : output real : Target gain
!! etadrv : output real : Driver efficiency
!! This routine calculates the parameters of a heavy ion driver
!! suitable for inertial fusion energy.
!! <P>Currently, the complicated model taken from the reference is not
!! complete, so it is recommended that the simple model is used
!! (set <CODE>ISIMP=1</CODE>)
!! Heavy-ion Driver Design and Scaling, R. Bieri et al.,
!! Fusion Technology, vol.21 (1992) 1583
!! Meier and Bieri, Fusion Technology, vol.21 (1992) 1547
!!
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use constants, only: echarge, mproton, degrad, pi
implicit none
! Arguments
real(dp), intent(in) :: aaion,bmax,dpp,dtheta,edrive,emitt,etai,lf,qion, &
sigma,sigma0,tauf,theta,vi
integer, intent(in) :: nbeams
real(dp), intent(out) :: gain,etadrv
! Local variables
real(dp), parameter :: c2 = 8.98755178737D16
integer, parameter :: isimp = 1 ! Switch for simple model (1=yes)
real(dp) :: ci,de,dgap,dlcore,drcore,e,eomc2,fins,floss, &
ibf,ibfo,ibi,ibpc,lfocus,lpf,lpfo,lpi,lppc,lq,phif,phifo, &
phii,phipc,rion,rs,rs1,rs2,rs3,rs4,sig,sig0,taufo,taui, &
taupc,tbrad,vf,vfo,vpc,vs,vscore,vshe,vsle,vspc,xhe,xle,xpc
real(dp), dimension(10) :: eve,gve
integer :: i,ie,j,nche,ncle,ncpc,nqarrs,nqle,nqpche,nquads
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (isimp /= 1) then
! Complex IONDRV model (incomplete)...
! ====================
! Electron charge / (proton mass * c**2)
eomc2 = echarge / (mproton*c2)
! Degrees to radians
sig = sigma * degrad
sig0 = sigma0 * degrad
! Beam current coefficient
ci = cbeam(etai,bmax,aaion,qion,sig,sig0,emitt,vi,eomc2)
! Final voltage (V)
vf = (edrive/(dble(nbeams)*tauf*ci))**(0.666666D0)
! Final beam current (A)
ibf = ci*sqrt(vf)
! Ion range (g/cm^2)
rion = (3.04D-5 * aaion + 349.0D0*aaion**(-2.21D0)) * &
(qion*vf/1.0D9)**(1.72D0-0.00275D0*aaion)
! Convert to kg/m^2
rion = rion*10.0D0
! Total spot size (rms of four terms):
! (1) Emittance term
rs1 = emitt / (betgam(aaion,qion,vf)*theta)
! (2) Dispersion (momentum spread) term
rs2 = 8.0D0*dpp*lf*theta
! (3) Alignment and jitter term
rs3 = dtheta*lf
! (4) Space charge term (inaccurate at present)
rs4 = 1.0D-3
rs = sqrt(rs1**2 + rs2**2 + rs3**2 + rs4**2)
! Target gain (inaccurate at present)
gain = 62.6D0 - 142.3D0*rs**1.25D0 * rion**0.633D0 + &
(32.8D0+82.1D0*rs)**0.783D0 * rion**0.647D0 * &
log(1.0D-6*edrive)
! Source/Injector characteristics:
! Transport beam radius (m)
tbrad = ( (emitt**2 * sig0 * c2)/(sig**2 * bmax/1.5D0 * etai * &
sqrt(vi)) )**(0.333333D0) * sqrt(2.0D0 * aaion * mproton / &
(qion * echarge))
! Extractor voltage (V)
vs = 4.0D0/3.0D0 * vi
! Source extraction gap width (m)
dgap = 0.8D0
! Injected current per beam (A)
ibi = 4.0D0 * pi * tbrad**2 * 5.46D-8 * &
sqrt( qion*vs**3/(aaion*dgap**4) )
! Beam loss fraction between injector and target
floss = 0.08D0
! Pulse duration at injection (s)
taui = 1.0D0/(1.0D0-floss) * tauf * ibf / ibi
! Pulse length at injection (m)
lpi = taui * sqrt( 2.0D0 * qion * echarge * vi / (aaion * mproton) )
! Initial voltage gradient (V/m)
phii = 0.3D0 * vi / lpi
! End of Low Energy Transport Stage:
! Insulated fraction of driver length
fins = 0.85D0
! Maximum voltage gradient, limited by insulator flash-over (V/m)
phifo = 1.0D6 * fins
! Voltage at which this limit is reached (V)
! The end of the low energy stage is set by the point at which
! the limit is reached.
vfo = (phifo/phii)**(2.0D0/3.0D0) * vi
! Beam current (A)
ibfo = ibi * (vfo/vi)
! Pulse duration (s)
taufo = taui * (vi/vfo)
! Pulse length (m)
lpfo = taui * vi / vfo * sqrt( 2.0D0 * qion * echarge * vfo / &
(aaion * mproton) )
! Length of the low energy transport stage (m)
! (rearrangement of integral of PHI(V).dl)
xle = 2.0D0*vi**1.5D0 / phii * (1.0D0/sqrt(vi) - 1.0D0/sqrt(vfo))
! End of Pulse Compression Stage:
! Pulse length = final pulse length (m)
lpf = tauf * sqrt( 2.0D0 * qion * echarge * vf / (aaion * mproton) )
lppc = lpf
! Length of the pulse compression region (m)
xpc = (lpfo - lpf)/0.3D0
! Voltage gradient (V/m) - limit already reached
phipc = phifo
! Voltage (V)
vpc = vfo + (xpc*phipc)
! Beam current (A)
ibpc = ibfo * (lpfo/lppc) * sqrt(vpc/vfo)
! Pulse duration (s)
taupc = taufo * (lppc/lpfo) * sqrt(vfo/vpc)
! End of High Energy Transport Stage:
! Voltage gradient (V/m) - limit already reached
phif = phifo
! Length of the high energy transport stage (m)
xhe = (vf - vfo)/phifo
write(*,*) ' Injector:'
write(*,*) ' vi = ',vi
write(*,*) ' ibi = ',ibi
write(*,*) ' taui = ',taui
write(*,*) ' phii = ',phii
write(*,*) ' lpi = ',lpi
write(*,*) ' End of LET stage:'
write(*,*) ' vfo = ',vfo
write(*,*) ' ibfo = ',ibfo
write(*,*) 'taufo = ',taufo
write(*,*) 'phifo = ',phifo
write(*,*) ' lpfo = ',lpfo
write(*,*) ' xle = ',xle
write(*,*) ' End of PC stage:'
write(*,*) ' vpc = ',vpc
write(*,*) ' ibpc = ',ibpc
write(*,*) 'taupc = ',taupc
write(*,*) 'phipc = ',phipc
write(*,*) ' lppc = ',lppc
write(*,*) ' xpc = ',xpc
write(*,*) ' End of HET stage:'
write(*,*) ' vf = ',vf
write(*,*) ' ibf = ',ibf
write(*,*) ' tauf = ',tauf
write(*,*) ' phif = ',phif
write(*,*) ' lpf = ',lpf
write(*,*) ' xhe = ',xhe
! Volt-second requirements and number of cores
! Volt-seconds per core (2.4T is field swing of the core
! material metglas)
drcore = 0.8D0
dlcore = 0.2D0
vscore = 2.4D0 * drcore * dlcore
! Low energy stage
vsle = taui * vi * log(vfo/vi)
ncle = 1 + int(vsle/vscore)
! Pulse compression stage
vspc = 2.0D0 * tauf * sqrt(vf) * (sqrt(vpc)-sqrt(vfo))
ncpc = 1 + int(vspc/vscore)
! High energy stage
drcore = 0.4D0
dlcore = 0.1D0
vscore = 2.4D0 * DRCORE * DLCORE
vshe = 2.0D0 * tauf * vf * (1.0D0 - sqrt(vpc/vf))
nche = 1 + int(vshe/vscore)
write(*,*) ' Number of cores: '
write(*,*) ' ncle = ',ncle
write(*,*) ' ncpc = ',ncpc
write(*,*) ' nche = ',nche
! Quadrupole requirements
! Effective quadrupole length (m)
lq = ( emitt * etai * sig0**2 * sqrt(c2*vi) / &
(sig*(bmax/1.5D0)**2) )**(0.333333D0) * &
sqrt(2.0D0*aaion*mproton/(qion*echarge))
write(*,*) ' lq = ',lq
! Actual quadrupole length (m)
write(*,*) ' alq = ',lq/etai
! Length of quadrupole focussing fields: low energy stage
lfocus = etai * vi/phii * (1.0D0 - vi/vfo)
nqle = 1 + int(lfocus/lq)
! Pulse compression + high energy stages
lfocus = 2.0D0 * etai * sqrt(vi)/phifo * (sqrt(vf) - sqrt(vfo))
nqpche = 1 + int(lfocus/lq)
! Number of quadrupole arrays
nqarrs = nqle + nqpche
! Number of quadrupoles
nquads = nqarrs * nbeams
write(*,*) 'nqarrs = ',nqarrs
write(*,*) 'nquads = ',nquads
etadrv = 0.28D0
else
! Simple model
! ============
! gve(k): target gain at edrive = k MegaJoules
gve(1) = 25.0D0
gve(2) = 44.0D0
gve(3) = 62.0D0
gve(4) = 76.0D0
gve(5) = 87.0D0
gve(6) = 97.0D0
gve(7) = 107.0D0
gve(8) = 115.0D0
gve(9) = 125.0D0
gve(10) = 132.0D0
! eve(k): driver efficiency at edrive = k MegaJoules
eve(1) = 0.232D0
eve(2) = 0.256D0
eve(3) = 0.269D0
eve(4) = 0.276D0
eve(5) = 0.282D0
eve(6) = 0.286D0
eve(7) = 0.290D0
eve(8) = 0.292D0
eve(9) = 0.294D0
eve(10) = 0.296D0
e = 1.0D-6 * edrive
ie = int(e)
de = e - dble(IE)
! Assume linear interpolations and extrapolations
! Would be better to prevent extrapolation
if (ie <= 1) then
gain = gve(2) - 1.0D-6*(edrive-2.0D6)*(gve(1)-gve(2))
etadrv = eve(2) - 1.0D-6*(edrive-2.0d6)*(eve(1)-eve(2))
else if (ie >= 9) then
gain = gve(9) + 1.0D-6*(edrive-9.0D6)*(gve(10)-gve(9))
etadrv = eve(9) + 1.0D-6*(edrive-9.0D6)*(eve(10)-eve(9))
else
gain = gve(ie) + de*(gve(ie+1)-gve(ie))
etadrv = eve(ie) + de*(eve(ie+1)-eve(ie))
end if
! Ensure sensible values
gain = max(0.01D0,gain)
etadrv = max(0.01D0,etadrv)
end if
contains
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function cbeam(etai,bmax,aaion,qion,sigma,sigma0,emitt,vi,eomc2)
!! Routine to evaluate the beam current coefficient, CI
!! author: P J Knight, CCFE, Culham Science Centre
!! etai : input real : Axial quadrupole packing fraction at injection
!! bmax : input real : Maximum field at the superconductor (T)
!! aaion : input real : Ion mass (amu)
!! qion : input real : Ion charge state
!! sigma : input real : Depressed tune (incl. space charge effects) (rad)
!! sigma0 : input real : Phase advance per lattice period (tune) (rad)
!! emitt : input real : Normalised beam emittance (metre-rad)
!! vi : input real : Injection voltage (V)
!! eomc2 : input real : Electron charge / (proton mass * light speed**2)
!! This routine calculates the beam current coefficient
!! of a heavy ion driver suitable for inertial fusion energy.
!! Heavy-ion Driver Design and Scaling, R. Bieri et al.,
!! Fusion Technology, vol.21 (1992) 1583
!! !
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
real(dp) :: cbeam
! Arguments
real(dp), intent(in) :: etai,bmax,aaion,qion,sigma,sigma0, &
emitt,vi,eomc2
! Local variables
real(dp) :: bbe
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Maximum field at the beam edge
bbe = bmax/1.5D0
! Perform calculation
cbeam = 2.89D6*(1.0D0-(sigma/sigma0)**2) * &
( sigma0**4 * etai**2 * vi * aaion/qion * (emitt/sigma)**2 &
* bbe**2 )**0.333333D0 &
* (2.0D0 * qion * eomc2 / aaion)**(5.0D0/6.0D0)
end function cbeam
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function betgam(aaion,qion,v)
!! Routine to calculate the relativistic factor (beta*gamma) for the
!! heavy ions in the beam driver
!! author: P J Knight, CCFE, Culham Science Centre
!! aaion : input real : Ion mass (amu)
!! qion : input real : Ion charge state
!! v : input real : Acceleration voltage (V)
!! This routine calculates the relativistic factor (beta*gamma) for the
!! heavy ions in a beam driver suitable for inertial fusion energy.
!! Heavy-ion Driver Design and Scaling, R. Bieri et al.,
!! Fusion Technology, vol.21 (1992) 1583
!! !
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
real(dp) :: betgam
! Arguments
real(dp), intent(in) :: aaion,qion,v
! Local variables
real(dp), parameter :: c2 = 8.98755178737D16
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
betgam = sqrt( 2.0D0*qion*echarge*v / (aaion*mproton*c2) )
end function betgam
end subroutine iondrv