fusion-product-sampling/reactivities.py at master · ukaea/fusion-product-sampling · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
from math import exp
import numpy as np

def bosch_hale_reactivity( Ti ):
    '''
    Evaluates the Bosch-Hale reactivity. Valid from 0.2 - 100keV with 0.25% accuracy.
    Note the BH fit returns the reactivity in cm**3 / s hence the conversion at the end.
    https://iopscience.iop.org/article/10.1088/0029-5515/32/4/I07/pdf
    '''

    # Coefficients
    C1 = 1.17302e-9
    C2 = 1.51361e-2
    C3 = 7.51886e-2
    C4 = 4.60643e-3
    C5 = 1.35e-2
    C6 =-1.06750e-4
    C7 = 1.366e-5

    # Gamov factor and reduced mass * c**2
    BG   = 34.3827   # [sqrt[keV]]
    mrc2 = 1124656   # [keV]

    # Get Theta
    theta = Ti / ( 1.0 - ( Ti * ( C2 + Ti * ( C4 + Ti * C6 ) ) ) / ( 1.0 + Ti * ( C3 + Ti * ( C5 + Ti * C7 ) ) ) )

    # Get Chi
    chi = ( BG**2.0 / (4.0 * theta) )**(1.0/3.0)

    # Get reactivity
    R = C1 * theta * ( chi / ( mrc2 * Ti**3.0 ) )**0.5
    R = R * exp( -3.0 * chi )

    # Convert to m**3 / s
    R = R / 1.0e6

    return R

def sadler_van_belle_reactivity( Ti ):
    '''
    Evaluates the Sadler-Van-Belle formula for the reactivity - Ti is temperature in keV.
    '''

    C1 = 2.5663271e-18
    C2 = 19.983026
    C3 = 2.5077133e-2
    C4 = 2.5773408e-3
    C5 = 6.1880463e-5
    C6 = 6.6024089e-2
    C7 = 8.1215505e-3

    U = 1.0 - Ti * ( C3 + Ti * (C4 - C5*Ti)) / ( 1.0 + Ti * (C6 + C7*Ti))

    R = ( C1 / ( U**(5.0/6.0) * Ti**(2.0/3.0) ) ) * exp( -1.0 * C2 * ( U / Ti )**(1.0/3.0))
    return R

def rate( ni, Ti, params ):
    '''
    Evaluates volumetric reaction rate.
    '''

    if( params['bosch-hale'] ):
        R = ni**2.0 * bosch_hale_reactivity(Ti)
    else:
        R = ni**2.0 * sadler_van_belle_reactivity(Ti)

    return R