УДК 539.173.84
АЛГОРИТМ МОДЕЛИРОВАНИЯ ПРОЦЕССА ПРОХОЖДЕНИЯ НЕЙТРОНОВ ЧЕРЕЗ ВЕЩЕСТВО
ПРОХОРЕЦИ.М., ПРОХОРЕЦ С.И., ХАЖМУРАДОВ М.А.
Описывается алгоритм расчета прохождения нейтронов через среду. Приводятся примеры расчета прохождения моноэнергетического пучка нейтронов через плоский и цилиндрический слои полиэтилена.
Введение
Последние годы отмечены расширением сферы использования нейтронов для решения различных научных и прикладных задач, ростом потребностей в источниках нейтронов и развитием математических методов описания взаимодействия этих частиц со средой [1]. Среди глобальных задач, в которых возникает необходимость учета взаимодействия нейтронов с веществом, можно отметить исследование и разработку ядерных энергетических систем с естественной безопасностью [2], сжигание отходов ядерного топлива быстрыми нейтронами [3], экологические проблемы после аварии на ЧАЭС [4]. Математические программы и алгоритмы, описывающие поведение нейтронов в таких сложных системах, не являются по многим причинам общедоступным интеллектуальным продуктом в Украине.
1. Цель работы
Целью настоящего исследования является разработка на основе ранее рассмотренной авторами модели взаимодействия нейтронов с веществом [5] алгоритма для расчета прохождения нейтронов через тела различной геометрической формы и различного химического состава. На практике для численной реализации таких задач используется статистический метод Монте-Карло.
2. Алгоритм прохождения нейтронов через однородную среду
В качестве примера рассмотрим задачу прохождения нейтронов через область f(x,y,z), заполненную однородным веществом (например, полиэтилен С2Н 4). Приведем основные этапы моделирования (рис.1):
1. Определяются геометрические размеры области f(x,y,z) в “мировой” системе координат X, Y,Z :
xmax _ a ; y max _ b ; zmax _ c ; Пр°изв°дИТся ввод и
подготовка других исходных данных.
2. Блок 2 представляет собой начало цикла по заданному числу разыгрываемых траекторий нейтронов.
3. В блоке 3 задается фиксированное значение начальной энергии нейтрона Б0 или разыгрывается энергия нейтрона б в падающем потоке.
РИ, 2003, № 4
4. Определяются или разыгрываются начальные координаты и направления влета каждого нового нейтрона (блок 4).
5. Анализируется энергия нейтрона (блок 5). Если энергия нейтрона равна или меньше заданной тепловой энергии т, то все процессы разыгрываются по специальному алгоритму (блоки 25-27).
6. В блоке 6 производится расчет длины свободного пробега и координат взаимодействия нейтронов в области X, Y, Z . Как показано в [5], длина свободного пробега
х = —у, (1)
где случайное число 0 <у<1, а £ t — полное макроскопическое поперечное сечение рассеяния нейтрона, измеряемое в см-1 [6]. Чтобы получить значение макроскопического поперечного сечения, величину микроскопического сечения [5] необходимо умножить на плотность вещества. По определению ^ t равно суммарной вероятности взаимодействия с различными веществами замедлителя, поэтому
Z t =Е 1 +Z 2 + - + Z m =Еt(H) +Zt(C) , (2)
Здесь ^t (H) и S t (C) — полные макроскопические сечения на водороде и углероде.
7. Производится проверка на вылет нейтрона в окружающую среду (блок 7).
8. Выбирается ядро (блок 8), с которым нейтрон вступил во взаимодействие при условии, что суммарная вероятность взаимодействия с ядрами водорода и углерода равна единице. Тогда вероятность того, что нейтрон провзаимодействует с водородом
pH и углеродом рс , можно представить соотношениями
ph
Z t
I t(H)+E t(Cr
(3)
Рс
z t(C)
I t(H)+E t(C)
(4)
а рс + рн = 1. Если 0 < у < рн , то взаимодействие произойдет с водородом, если рн < у < рн + рс , то с ядром углерода, где у — очередное случайное число из датчика случайных чисел.
9.Разыгрывается тип взаимодействия нейтрона с выбранным ядром. Для ядра атома водорода таким взаимодействием будет упругое рассеяние или радиационный захват (блоки 11 и 12 соответственно). Так как
Z t(H) = Ze(H) + Е cap(H) (5)
и ре + рсар = 1, то вероятность упругого рассеяния
125
Рис.1. Блок - схема алгоритма программы расчета прохождения нейтронов через вещество: САЗ — счетчик актов захвата; САД — счетчик актов деления; УР — упругое рассеяние; КБ — кинематический блок; НУР — неупругое рассеяние; САВ — счетчик актов вылета
Ре = Е е/Е t (6) Если очередное случайное число 0 < у < 1 удовлетворяет условию 0 < у < ре , то будет упругое рассе-
и захвата рсар _ Е сар /Еt • (7) яние, а в случае Ре < У ^ Рсар + Ре имеет место захват нейтрона.
126
РИ, 2003, № 4
10. В блоке 13 происходит розыгрыш и определение кинетических параметров рассеянного нейтрона — углов и энергии, полагая А = 1 в формуле [5]. Используя очередное случайное число у , азимутальный угол можно представить в виде
Ф = 2лу1. (8)
При розыгрыше полярного угла в системе центра масс ядро - налетающий нейтрон воспользуемся тем, что если энергия нейтрона <0,1МэВ, угловое распределение нейтронов абсолютно изотропно по cos 9cm, а для энергий нейтронов больше 0,1МэВ можно воспользоваться таблицей [7], в которой приведено угловое распределение упруго рассеянных нейтронов р(р,Е) (р = cos9cm). При фиксированной энергии е вероятность р нормирована следующим образом:
с нейтроном. Поэтому этапы моделирования в этом случае будут происходить в такой последовательности.
— Разыгрываем тип взаимодействия: упругое и неупругое рассеяние, захват нейтрона (реакция n, у ) и деление ядра.
— При упругом или неупругом рассеянии разыгрываем новое направление нейтрона, возникшего в результате взаимодействия, и его энергию. Обычно предполагают, что неупругое рассеяние изотропно в лабораторной системе координат, поэтому для розыгрыша азимутального угла Ф используется выражение (8), а полярного угла в лабораторной системе — выражение
9 = 2у-1, (11)
где у — случайное число в интервале (0,1).
m=M
Е Pm (Pm) APm = 1 , (9)
m=1
где М — количество равновероятных групп по косинусу угла рассеяния 9cm, приведенных в таблице [7]. Поэтому для розыгрыша р извлекаем из генератора случайных чисел очередное 0 < у2 < 1 и, решив неравенство
р1ДР1 + р2ДР2 + ••• + pm-1APm-1 <72 <
< Р1ДР1 +12 + р2ДР2 + ••• pmДРш , (10)
находим р = р m. Осуществив переход в лабораторную систему координат, находим новую энергию в очередном акте рассеяния нейтрона [5].
11. По аналогии с рассеянием нейтрона на ядре водорода определяем тип взаимодействия с ядром углерода. В интервале с энергией налетающего нейтрона до 14МэВ полное сечение взаимодействия с углеродом определяется упругим рассеянием. В [7] угловое распределение упруго рассеянных нейтронов представлено в виде таблицы с шагом по косинусу угла рассеяния 0,1 и шагом по энергии 0,5МэВ.
12 .С использованием генератора случайных чисел для рассеянного на ядре углерода нейтрона находим азимутальный угол, косинус полярного угла в системе центра масс налетающий нейтрон - ядро 12 С и, наконец, энергию (по аналогии с п.10). Ввиду отсутствия у авторов данных для углерода в области энергий до 0,5МэВ при определении кинематических параметров предполагаем, что угловое распределение рассеянных нейтронов изотропно в системе центра масс.
3. Алгоритм прохождения нейтронов через среду с делящимся материалом
При наличии в среде делящегося материала, например 235 и , в рассмотренном выше алгоритме появляется дополнительная ветвь вычислений (блоки 20 — 24, 30 — 33), которая начинается с выбора ядра делящегося материала как объекта взаимодействия
—Энергию каждого неупругого рассеянного нейтрона также разыгрываем. Если распределение энергии е' рассеянных нейтронов можно считать максвелловским, то формулу для розыгрышае' приводим к виду [6]
1 -ГЕ'/T(1 + E'/T) = у 1 -Гь/1(1 + E/T) , (12)
,-E/T/
здесь T = aVE (a — постоянная, зависящая от рода вещества); е — энергия налетающего нейтрона, а
0 < Е' < Е; у — случайное число в диапазоне (0,1).
— Сравниваем энергию выходящего нейтрона с заданной во вводном блоке пороговой энергией
En , равной энергии связи нейтрона в ядре. Если
Е' > Еn, то переходим к п.2 алгоритма. Если Е' < En, то неупругие процессы не разыгрываются. Если нейтрон претерпел радиационный захват, то расчет траектории нейтрона завершается. Если число разыгранных траекторий меньше заданного, то переходим к п.2 алгоритма.
— Разыгрываем число возникших в реакции деления нейтронов v . Это случайная величина, для которой известно ее среднее значение V . Для ее розыгрыша вводим дополнительное целое число m и полагаем m < v < m +1. Также полагаем, что v может принимать два значения: V1 = m и V2 = m +1 с вероятностями P1 = m +1 -v и P2 =v-m [6]. Следовательно, если случайное число v < P1, то необходимо считать v = V1, а если P1 < v, то считать V = V2 .
— Из отношения количества запаздывающих vd к полному (суммарному) количеству нейтронов vt определяем тип возникшего нейтрона — мгновенный или запаздывающий. Если
y<vd/ Vt, (13)
то возникает запаздывающий нейтрон, а если
РИ, 2003, № 4
127
y> vd/ vt
(14)
то появляется мгновенный нейтрон.
—Разыгрываем углы вылета всех мгновенных нейтронов, полагая, что для каждого нейтрона деления все направления вылета равновероятны. Определяем энергию каждого мгновенного нейтрона, используя соотношения [5].
—Разыгрываем энергию, углы и время вылета запаздывающих нейтронов. Если в возникшей ветви дерева энергия быстрого или запаздывающего нейтрона меньше En , то неупругие процессы (и деление в том числе) не рассматриваются. Если Е' > En , то переходим к п.2 алгоритма.
Таким образом, для расчета взаимодействия нейтронов с энергией до 14МэВ со средой необходимо иметь целый набор экспериментальных данных для большого числа химических элементов. В это число входят данные по полным сечениям упругого и неупругого рассеяния, поглощения, радиационного захвата, ядерных реакций и т.п. Из этого перечня видно, что потребность в ядерных данных для расчета прохождения нейтронов через полиэтилен может удовлетворить справочник, в котором приведены константы взаимодействия нейтронов с элементами, входящими в состав атмосферы и земной коры [7]. Однако для проведения более сложных и точных расчетов, в частности, при наличии делящихся материалов необходимо пользоваться библиотеками оцененных ядерных данных ENDF/B-IV или ENDL.
4.Примеры численной реализации и сравнение с литературными данными
Наличие данных о нейтронных сечениях для ядер углерода и водорода дает возможность провести численный расчет прохождения нейтронов через плоский и цилиндрический слои полиэтилена. Решение этой задачи позволяет проверить работу алгоритма программы путем сравнения расчетных величин с имеющимися литературными данными
[8,9].
На рис. 2 и 3 представлены спектр выходящих из
полиэтиленовой пластины нейтронов AN/ Ди и распределение источников тепловых нейтронов
S(Z), как функция толщины слоя пластины Z .и —летаргия нейтронов, эта величина равна среднему числу соударений, которое необходимо для замедления нейтронов от энергии Е0 = 2МэВ до тепловой энергии. Источниками тепловых нейтронов S(Z) служат последние их столкновения перед тем, как они превратились в тепловые.
Сравнение результатов расчета с использованием разработанного алгоритма и имеющимися литературными данными [8] для случая плоской пластины полиэтилена указывает на их согласие в пределах около + 20%.
I!
Рис. 2.Спектр нейтронов с начальной энергией Е0 = 2МэВ на выходе из пластины полиэтилена толщиной 5см [8]
S
140
120
100
80
60
40
Рис. З.Распределение источников тепловых нейтронов (Е < 107 МэВ ) по толщине пластины (Н = 5см, Е0 = 2МэВ) [8]
Методом Монте-Карло также было рассчитано прохождение нейтронов с начальными энергиями до ЗМэВ через цилиндрические слои полиэтилена с внешними диаметрами до 120см. Расчет проводился для случая падения узкого пучка нейтронов, направленного нормально к оси цилиндра по его диаметру. Энергетическое распределение нейтронов, прошедших во внутреннюю область цилиндра, приведено на рис.4. Полученные нами результаты расчета для цилиндра с внешним диаметром 120см согласуются с данными [9]. Они подтверждают также вывод работы [9] о зависимости ослабления потока нейтронов не только от толщины, но и от диаметра цилиндрического слоя. Следовательно, этот фактор необходимо учитывать при проектировании спектрометров и дозиметров нейтронов.
128
РИ, 2003, № 4
Рис.4. Спектр нейтронов с начальной энергией Eg = 3МэВ внутри цилиндра с внешним диаметром: а — 120см (толщина стенки 9см); б — 5,5см (толщина стенки 1 см)
Заключение
Полученные в данной работе результаты свидетельствуют о работоспособности алгоритма программы для расчета прохождения быстрых нейтронов через объекты различной геометрической формы и различного состава вещества. Он может быть принят за основу при создании математических программ для расчета сложных по форме и химическому составу объектов исследования и практического использования для нейтронной радиографии, экологически безопасной ядерной энергетики будущего, дозиметрии и т.п.
Литература: 1. Batiy V. G, ProkhoretsI.M., Prokhorets S.I., Khazhmuradov M.A. et al. Development of mathematical and experimental model of neutron radiography set up // Problems of Atomic Science and Technology. Series: Nuclear Physics Investigations. 2003. N2. P. 116-118. 2. Бойко B.A., Карнаухов И.М., Лапшин В.И. Усилитель мощности — основа ядерной энергии XXI века. Препринт ХФТИ. Харьков: ННЦ ХФТИ. 2001. 50с. 3. Бойко
B. A., Егоров A.M., Зайцев Б.В., Карнаухов И.М., Кобец A. Ф, Лапшин В.И. Сжигание отходов ядерного топлива быстрыми нейтронами в электроядерной энергетической установке — альтернатива геологическому захоронению: случай Украины. Препринт ХФТИ-2001-2. Харьков: ННЦ ХФТИ. 2001-2. 34с. 4. Слабоспицкий Р.П, Кочетов С.С., Кириченко В.В. Система экспрессной сортировки РАО с использованием ЛУЭ (линейного ускорителя электронов). Препринт ХФТИ 20021. Харьков: ННЦ ХФТИ. 2002-1. 19с. 5. Прохорец И.М, Прохорец С.И., Хажмурадов M.A. Математические модели прохождения нейтронов через вещество / / Радиоэлектроника и информатика. 2003. №1. С. 124-128 6. БусленкоН.П., Голенко Д.И, СобольИ.М., СраговичВ.Г., Шрейдер ЮЛ. Метод статических испытаний. М.: Физ-матгиз, 1962. 321с. 7. Медведев Ю.A., Степанов Б.М., Труханов Г.Я. Ядерно-физические константы взаимодействия нейтронов с элементами, входящими в состав атмосферы и земной коры. Справочник. М.: Энер-гоиздат, 1981. 303с. 8. Ермаков С.М., Золотухин В.Г., Петров Э.Е. Расчет прохождения нейтронов через плоский слой полиэтилена // Атомная энергия. 1963. Т. 15. Вып. 3. С. 253-255. 9. Мирошников Г.В., Мирошни-кова A.H., Солодухин H.A. Расчет прохождения нейтронов через цилиндрические слои полиэтилена методом Монте-Карло // Атомная энергия. 1969. Т.26. Вып.4.
C. 372-373.
Поступила в редколлегию 11.04.2003
Рецензент: д-р техн. наук, проф. Кривуля Г.Ф.
Прохорец Иван Михайлович, канд. физ.-мат. наук, старший научный сотрудник ННЦ ХФТИ. Научные интересы: информационно-измерительные системы и детекторы в ядерной физике. Адрес: Украина, 61108, Харьков, ул. Академическая, 1, тел. 35-61-13.
Прохорец Светлана Ивановна, инженер-исследователь Национального ННЦ ХФТИ. Научные интересы: математическое моделирование физических процессов и систем, программирование. Адрес: Украина, 61108, Харьков, ул. Академическая, 1, тел. 35-65-94.
Хажмурадов Манап Ахмадович, д-р техн. наук, профессор, начальник отдела ННЦ ХФТИ. Научные интересы: математическое моделирование физических процессов и систем, автоматизация проектирования, программирование. Адрес: Украина, 61108, Харьков, ул. Академическая, 1, тел. 35-68-46.
РИ, 2003, № 4
129