CC BY
53
ФИЗИКА АТОМНОГО ЯДРА И ЭЛЕМЕНТАРНЫХ ЧАСТИЦ
Математическая модель рассеяния фотонов в веществе в задачах расчета и оптимизации радиационной защиты для инспекционно-досмотровых комплексов
В. И. Петрунинa, С. А. Огородниковb, М. А. Арлычев, И. Е. Шевелев
ООО Лаборатория Скантроник. Россия, 190013, Санкт-Петербург, ул. Введенского канала, д. 7.
E-mail: a [email protected], b [email protected] Статья поступила 17.11.2014, подписана в печать 02.12.2014.
Излагается математический аппарат оригинального программного обеспечения, разработанного авторами для расчета радиационной защиты и полей тормозного и рассеянного излучения линейных ускорителей электронов, которое адаптировано к широкому классу прикладных задач. В основе программы лежит метод Монте-Карло, позволяющий проводить статистический расчет потока квантов тормозного излучения и выбирать траектории квантов тормозного излучения при помощи генератора случайных чисел. Программа предназначена для проектирования систем формирования пучков и полей излучения, защит сложного профиля с лабиринтами и каналами, оптимизации размеров защитных барьеров, для оценки ожидаемого отклика детекторной системы и распределения дозы по сечению объекта. Для верификации достоверности расчетов, проводимых с помощью оригинального программного обеспечения, в статье приводятся сравнения результатов численных экспериментов с помощью кода, предложенного авторами, и Монте-Карло кода MCNP версии 4C3 для обратно рассеянных и прямо прошедших квантов тормозного излучения с граничной энергией 24 МэВ для ряда материалов и толщин.
Ключевые слова: тормозное излучение, поле излучения, радиационная защита, ускоритель электронов, Монте-Карло, инспекционно-досмотровый комплекс, рассеяние, траектория, фотон, доза, прозрачность, альбедо.
УДК: 539.1.04. PACS: 28.41.Qb.
Введение
В мире создано несколько программ статистического моделирования прохождения гамма-квантов через вещество, однако до сих пор не удалось создать достаточно универсального кода на «все случаи жизни». В крупных исследовательских лабораториях создаются собственные программы MCNP, GEANT4 [1, 2], есть также программы Penelope, Fluka, EGS [3-5]. Эти программы имеют специфику, связанную с тем, что они нацелены в первую очередь на решение определенных задач — расчет систем детектирования частиц в ядерной физике и физике высоких энергий, расчет поглощенных доз и планирование облучения в медицине, расчет защиты ядерных реакторов и т. п. Безусловно, эти программы могут быть использованы и для расчета радиационной защиты. Однако сложность конструкции инспекционно-досмотровых комплексов и специфика их размещения затрудняют использование перечисленных программ для расчета защиты этих комплексов.
Поэтому для решения широкого круга прикладных задач, связанных с расчетом и оптимизацией радиационной защиты инспекционно-досмотровых комплексов (ИДК) различных классов на основе ускорителей электронов и особенно для тех типов, в которых пучок тормозного излучения выводится в окружающее пространство (например, мобильный, портальный или стационарный железнодорожный ИДК), необходим не
только специализированный, но и достаточно точный и быстрый программный инструмент расчета.
Аналогичная ситуация сложилась при проектировании стационарной радиационной защиты для ускорительных комплексов для различных прикладных целей (например, стерилизаторов и т.п.), где ослабление излучения достигается многократным рассеянием квантов в защитных стенах тоннелей.
Кроме того, в различных задачах интроскопии необходим точный расчет и оптимизация трактов коллимации пучков высокоэнергетичных фотонов обеспечения высоких характеристик комплексов.
Все вышесказанное вынудило авторов провести собственную разработку кода, адаптированного к широкому классу прикладных задач.
В основе программы лежит метод Монте-Карло, позволяющий проводить статистический расчет потока квантов тормозного излучения и выбирать траектории квантов тормозного излучения при помощи генератора случайных чисел.
В специальном программном обеспечении учитывается как математический формульный аппарат процессов генерации тормозного фотонного излучения в веществе, так и табулированные данные массового коэффициента ослабления гамма квантов в веществе.
Программа предназначена для проектирования систем формирования пучков и полей излучения, защит сложного профиля с лабиринтами и каналами, опти-
мизации размеров защитных барьеров. Кроме того, результаты расчета могут быть использованы для оценки ожидаемого отклика детекторной системы и распределения дозы по сечению объекта.
В программе предусмотрен метод «конструирования» барьеров на экране дисплея в двух координатах, например у—г в виде изображения растрового формата, с одновременной (не визуальной) установкой третьей координаты, причем гетерогенность барьера определяется различными оттенками цвета. Предусмотрены также визуализация прохождения квантов через барьер, построение графиков углового и энергетического распределения прямо прошедших, рассеянных и отраженных квантов.
Необходимо заметить, что излучение, генерируемое при торможении вторичных электронов и позитронов, образующихся вследствие фотоэффекта и эффекта рождения электрон-позитронных пар, не учитывалось, так как для энергетического диапазона тормозного излучения, который используется при расчете и проектировании инспекционно-досмотровых комплексов, и в котором превалирует комптон-эффект, вкладом данных эффектов можно пренебречь.
1. Математическая модель 1.1. Энергетическое и угловое распределение фотонов тормозного излучения
Энергетическое и угловое распределение фотонов тормозного излучения, генерируемого в мишенях ускорителей, зависит от ряда параметров — материала и толщины мишеней, спектрального и углового распределения падающих на мишень ускоренных электронов, материала и конструкции мишенного узла и первичного коллиматора. Определяющим в энергетическом и угловом распределении фотонов тормозного излучения является двойное дифференциальное по энергии и углу вылета сечение образования тормозного излучения, имеющее достаточно сложный вид. Обычно в полуаналитических методах расчета радиационной защиты используются модели с разделением функций распределения фотонов по энергиям и углу вылета, так как точный расчет двойного дифференциального по энергии и углу распределения квантов тормозного излучения является уточняющей поправкой более высокого порядка, чем аппроксимационные спектральные и угловые распределения.
Спектральное распределение потока энергии тормозного излучения определяется по формуле Шиффа [8]
dE. {Ee'E7'Z
= 8 • 2 • (1 - Krei) • (ln a - 1)+ Kr2el • ln a - 2 , (1)
(in«- 2))
где
a1 = 2 •
/
(2)
Ee + Me (1 - Krel)
Me
Krel
E
Krel = -EL
2 2 O" 1 0"2
a21 + a22 , a2 = 191 • exp ^-1 ln Z^,
Здесь Е1 — энергия рожденного кванта тормозного излучения, Ее — энергия электрона, Ме — масса покоя электрона, Z — атомный номер материала мишени.
В выражении (2) коэффициент для толстой мишени (по крайней мере оптимальной для выхода тормозного излучения вперед) равен 191. С таким значением коэффициента функция, заданная выражением (1), хорошо аппроксимирует спектральное распределение [9].
Для получения относительного дифференциального распределения числа фотонов в зависимости от их энергии необходимо значение интенсивности в относительных единицах
^1ге1(Ет)
г =
разделить на EY. Тогда
d1rel (ey ) =
dEY
dN (E~)
dEY
rel
Так как факт испускания фотона с энергией Еу под углом в является случайным, необходимо разыграть эти параметры с использованием случайных чисел. Энергия фотона определяется из соотношения
^ 1
dN E) dE.
dEY
E''=Ee dN (Ey) 0
dEY
dEY
Интегрирование в числителе ведется до того значения энергии в верхнем пределе интеграла, пока все выражение не будет равно АЕу, где АЕу — случайное число при розыгрыше энергии фотона.
Функция распределения фотонов по углам для толстой мишени из материала с высоким атомным номером Z и для энергии ускоренных электронов от 2 до 50 МэВ в плоскости сечения пучка, совпадающей с осью потока электронов, аппроксимируется формулой [10]
Р(в) 1
F (в) =
P(0)
1 + ( ^
1 + I 1.73
1.4
По сути F(в) = -Ny-, т. е. является вероятностью испускания фотонов в элемент телесного угла -У в направлении в. Поэтому
dN(в)= F(в), -У = F(в) sin в dв -ф,
где ф — азимутальный угол.
Учитывая азимутальную симметрию испускания фотона, можно написать
dN (в) = 2nF (в) sin в dв.
Отсюда следует, что для угла вылета в
А(в) =
j F(в) sin в dв
вт'т втах
j F(в) sin в dв
Ee + Me'
Углы втШ и втах определяются, например, коллимиру-ющей системой. Азимутальный угол ф разыгрывается из условия ф = 2пА.
a=
в
Определив таким образом энергию и направление вылета фотона, можно перейти к расчету траектории отдельного фотона (кванта) также методом случайных испытаний.
1.2. Расчет отдельного трека кванта тормозного излучения
Схему расчета отдельного трека кванта можно построить следующим образом.
Сначала находится пробег кванта с энергией Е7. Если /л(Е7г) — полный линейный коэффициент поглощения в единицах 1/см, то пробег кванта будет находиться в интервале от г до г + Дг с вероятностью
Д V(Еу, г) = /(Е^)в-^(Е->2)гДг.
Здесь в-Р(Е12)г — вероятность пройти расстояние г без столкновений и /(Е7, X)Дг — вероятность испытать взаимодействие на отрезке Дг (0 < г < то). Разобьем интервал 0 ^ 1, в котором меняется V(Е7, г), на п равных отрезков Д^ = 1/п, тогда каждому интервалу изменения V от до + Д^М^ можно сопоставить интервалы Дги, которые будут разными при равных Д V. Случайные числа Аг распределены равновероятно в интервале 0 ^ 1, поэтому при разбиении интервала случайных чисел на п равных частей можно написать
Д V = ДАкг,
1.3. Розыгрыш типа взаимодействия кванта тормозного излучения с веществом
Далее необходимо разыграть тип взаимодействия кванта тормозного излучения с веществом.
В данном расчете процессы взаимодействия кванта тормозного излучения, такие как фотоэффект и образование пар, учитывались только как факт поглощения, поэтому исходя из соотношения
ЛКШ + ^оЬ + Мра^ = 1
или
JKNT _ 1 Jphoto + Jpair
т.е. АсотрЬп = 1 - АаЬзогЬ. Если Асо^оп — случайное число, большее
чем /кнт/л , то происходит поглощение кванта, и он выбывает из рассмотрения. В противном случае происходит комптоновское рассеяние и можно перейти к его рассмотрению.
Далее разыгрывается угол рассеяния кванта при комптоновском процессе. Если Аш — случайное число, то
Ч ш0
¡чй^ М 4k sin ш du
АШ —
4п
т. е.
или
jdBdV М dk sin 0 0
Числитель данного выражения имеет вид [10]
AAkr _ j,(EY,Z)e-j(E->Z)rДг, Ar _ C - e-j(E'>Z)r.
1
Ош _ Г
2а2 (1 + а - а cos ш)2
C _ 1 при r ^ <ж , т. е.
4 + 10а - 8а2 + а3 - (4 + 16а2 + 2а3) cos ш0 +
Ar _ 1 - e-M(EYZ)r,
+ (6а + 10а2 + а3) cos2 ш0 - 2а2 cos3 ш0
+
или
1 - Аг _ e-IJ'(EYZ)r.
+
/а2 - 2а - 2\
V а3 )
В связи со случайным значением Ar в пределах где а _ E7/0.511, а знаменатель
ln(1 + а - 2 cos ш0)
0 ^ 1 число 1 - Ar также случайно:
A* _ e-j(EyZ)r,
2nr,2
/1 + а
V а2
^ - iln(1 + 2а)
1 + 2а а
+
или
in a:
= -/(Е7г )г.
Вышеизложенное справедливо для однородного материала X. В общем случае материал барьера может быть распределен в пространстве г, в котором распространяется квант, по определенному закону. В программе гамма квант «узнает» характер материала, и тогда 1п А* будет равен сумме /(Е7г(г; ))Дг, или в пределе
1 1 + 3а
+ H-ln(1 + 2а) -
2а
(1 + 2а)2
Энергия кванта после рассеяния определяется энергией налетающего кванта и углом рассеяния, определенном по соотношению
1 + а(1 - cos ш) '
rbegin
in a: _ -
j(Ey, r) dr.
геп(
Интегрирование ведется до верхнего предела геп(], пока логарифм случайного числа не станет равным интегралу, При розыгрыше 1-го пробега может оказаться, что геп( > гаЬкогЬ. Это означает, что квант прошел поглотитель без взаимодействия. В противном случае «история» кванта прослеживается дальше.
Здесь необходимо остановиться на минимальном значении энергии после рассеяния. Если а' окажется меньше установленного значения ат„, то история кванта заканчивается. Выбор ат„ требует некоторых замечаний. Для того чтобы проследить историю кванта с начальной энергией 1 МэВ до потери им энергии до 0.025 МэВ, требуется примерно в два раза больше рассеяний, чем до потери им энергии до 0.05 МэВ. Это существенно сказывается на скорости счета. В то же время установление высокого значения энергии «обрезания» («cut-off energy») может привести к ошибкам в расчете.
х
х
а
а_
Значение атш, которое не дает существенных ошибок, зависит от материала барьера. Рост коэффициента фотоэлектрического поглощения (а только этот процесс идет при малых энергиях) с увеличением атомного номера материала Z позволяет увеличить значение атш при больших Z. Численные эксперименты показали, что для водорода атш = 0.010 МэВ, для воды и бетона 0.025 МэВ, для железа 0.05 МэВ, для свинца 0.1 МэВ.
Азимутальный угол рассеяния = 2-кАх (см. угловые распределения, п. 1.1.)
После проведения операций по розыгрышу рассеяния можно вычислить состояние кванта после п-го рассеяния. Направление движения кванта в трехмерном случае в сферической системе координат записывается следующим образом:
СОЭ 6п+\ = С05б>„С05Ш„+1 +5т6>„5тш„+1 СОЭХп+Ь
cos(^„+i-^„) =
Sin0„+1 COS tjJtl+\ - COS 6tl COS 6tl+1 sin#„sin#„+1
Приращение координат для этих углов вычисляется по формулам
AZ„ = rn cos вп, АХп = rn sin вп cos (fn, AYn = rn sin вп sin (fn.
Таким образом, определено фазовое состояние кванта после и-го рассеяния.
2. Статистический расчет и сравнение с MCNP
Для верификации достоверности проводимых расчетов с помощью вышеизложенного математического аппарата в статье приводится сравнение расчета с помощью Монте-Карло кода MCNP версии 4СЗ [6] и данного специального программного обеспечения для обратно рассеянных и прямо прошедших квантов тормозного излучения с граничной энергией 24 МэВ в бетоне, железе и свинце для широкого диапазона толщин.
2.1. Прямо прошедшее излучение
Постановка расчетной геометрии состояла в задании бесконечно тонкого игольчатого пучка, падающего перпендикулярно поверхности экранирующего материала (рис. 1). В качестве экранирующих материалов использовались бетон (плотность р = 2.302 г/см3), железо (р = 7.874 г/см3) и свинец (р= 11.342 г/см3). Энергетический спектр тормозного излучения задавался апроксимационной формулой Шиффа (1), взятой в основу математической модели, описанной в настоящей статье. Толщины трех разных материалов варьировались с максимальными значениями для бетона 260 см, стали 60, и свинца 26 см. Для каждого эксперимента рассчитывалось 107 траекторий квантов тормозного излучения. Детектор брался протяженным (см. рис. 1), а энерговыделение в нем рассчитывалось в водяном эквиваленте.
Источник тормозного излучения (0,0,0)
Игольчатый пучок вдоль оси Z
Плоскость детектирования
Рис. 1. Геометрия для расчета прямо прошедшего излучения
2.2. Обратно рассеянное излучение
Геометрия для расчета обратно рассеянного излучения приведена на рис. 2. Падающий пучок брался игольчатой формы и падал перпендикулярно на поверхность материала. В моделировании подсчитывались обратно рассеянные от материала фотоны, пересекающие полусферу радиусом 1 м. Регистрирующий детектор определялся полярным углом с шагом 3.75°. В качестве рассеивающих материалов использовались бетон и железо толщинами 200 и 48 см соответственно. Альбедо рассчитывалось и нормировалось как дифференциальное дозовое альбедо с детектированием потока обратно рассеянных фотонов в водяном эквиваленте.
Падающий пучок
Обратно рассеяный пучок / .
/ п
Рис. 2. Геометрия для расчета обратно рассеянного излучения
2.3. Результаты расчетов
Результаты математического моделирования полной дозовой прозрачности барьеров из бетона, железа и свинца, сделанные с помощью рассматриваемой в данной статье математической модели, а также результаты расчетов с помощью МСЫР [6], приведены на рис. 3, 4 соответственно. В целом, несмотря на небольшое разногласие, результаты расчетов хорошо соотносятся друг с другом. Некоторое занижение расчетов согласно рассматриваемой математической модели, вероятно, связано с отсутствием учета вторичного тормозного излучения, а также аннигиляционного
1(Г
0 50 100 150 200 250 300 Толщина, см
0 10 20 30 40 50 60 Толщина, см
10 15 20 25 Толщина, см
30
Рис. 3. Полная дозовая прозрачность 24 МэВ тормозного излучения: а — в бетоне (плотность р = 2.302 г/см3); б - в железе (р = 7.874 г/см3); в — в свинце (р = 11.342 г/см3). Зеленым цветом — расчетные значения согласно рассматриваемой в статье математической модели. Синим цветом — расчетные значения с помощью MCNP [6]
излучения, связанного с рождением электрон-позитрон-ных пар в диапазоне энергий до 24 МэВ.
Результаты расчетов дифференциального дозового альбедо и сравнение с данными расчетов согласно MCNP [6] и данными IAEA Report [7] приведены на рис. 4. Из зависимостей, приведенных на графике, видно, что максимальное альбедо достигается при
о
ч
о ю л
10"-
s
I
о а
о
10-
Рис. 4. Дифференциальное дозовое альбедо для 24 МэВ тормозного излучения в виде игольчатого пучка, падающего нормально к плоскости рассеивающего материала. Данные включают в себя все фотоны, рассеянные в полярный угол ± 1.875°
строго обратном рассеянии излучения и заметно уменьшение альбедо с увеличением угла рассеяния ввиду того, что материал начинает больше экранировать выход фотонов. Кроме того, заметна разница: в результатах расчетов согласно MCNP дифференциальное дозовое альбедо для бетона и свинца практически совпадает, тогда как для расчетов согласно рассматриваемой в статье модели они довольно значительно разнятся, особенно под рассеянием строго назад.
Разногласие расчетов, возможно, связано с различным методом расчета и нормировки дифференциального дозового альбедо в математической модели, рассматриваемой в настоящей статье, а также, возможно, разным химическим составом бетона, взятого за основу в расчетах согласно MCNP и данной математической модели.
На наш взгляд, именно альбедо является наиболее чувствительным и показательным критерием для оценки математической модели и в целом видно неплохое совпадение результатов расчетов согласно MCNP [6], данных IAEA Report [7] и рассматриваемой модели, т. е. результаты расчетов согласно данной математической модели можно трактовать как достаточно достоверные.
Заключение
Референс расчеты для прямо прошедшего и обратно рассеянного излучения, сделанные с помощью специального программного обеспечения, основанного на математическом аппарате, изложенном в данной статье, показали хорошее согласие с результатами расчетов, сделанных с помощью расчетного кода MCNP версия 4СЗ [6], а также данными IAEA Report [7].
Математическая модель достоверно описывает процессы взаимодействия тормозного излучения с веще-
10"
0 10 20 30 40 50 60 70 80 90 Полярный угол испускания, град
— Бетон
- Бетон MCNP
■ Бетон IAEA
— Железо
- Железо MCNP
• Железо IAEA
ством и позволяет достаточно точно рассчитывать до-зовые поля нерассеянного и рассеянного излучения.
Расчетное ядро специального программного обеспечения и интерфейс для учета реальной трехмерной геометрии экранирующих и рассеивающих материалов позволяют проводить радиологический расчет систем с ускорителями электронов с условиями численного эксперимента, максимально приближенного к реальности.
Следует заметить, что специальное программное обеспечение является эффективным и удобным прикладным инструментом для оптимизации конфигурации радиационной защиты для инспекционно-до-смотровых комплексов и установок сканирования для различных целей на основе ускорителей электронов, для обеспечения радиационной безопасности в соответствии с требованиями российских действующих СанПиН 2.6.1.2369-08, ОСПОРБ-99/2010 и НРБ-99/2009.
Авторы глубоко признательны профессору, д.ф.-м.н. В. И. Шведунову за обсуждения, ценные замечания и комментарии.
Список литературы
1. https://mcnp.lanl.gov/
2. http://geant4.cern.ch/
3. http://www.oecd-nea.org/tools/abstract/detail/nea-1525
4. http://www.fluka.org/
5. http://www.nrc-cnrc.gc.ca/eng/solutions/advisory/ egsnrc_index.html
6. Avila-Rodriguez M.A., DeLuca P.M., Bohm T.D. // Radiation Protection Dosimetry. 2005. 116, N 1-4. P. 547.
7. International Atomic Energy Agency. Radiological safety aspects of the operation of electron linear accelerators. Rep. N 188. Vienna, IAEA, 1979.
8. Schiff L.I. // Phys. Rev. 1946. 70. P. 87.
9. Darby E.K. // Canad. J. Phys. 1951. 29(6). P. 569.
10. Стародубцев С.В., Романов А.Н. Взаимодействие гамма-излучения с веществом. Ч. 1. Ташкент, 1964.
11. Бурмистенко Ю.Н. Фотоядерный анализ состава вещества. М., 1986.
12. Петрунин В.И. Разработка систем таможенного и промышленного радиационного цифрового контроля крупногабаритных объектов на базе линейных электронных ускорителей: Дисс. ... докт. техн. наук. СПб., 2002.
A mathematical model of photon scattering in matter for problems of calculation and optimization of radiation shielding in X-ray inspection systems
V.I. Petrunina, S. A. Ogorodnikovb, M. A. Arlychev, I.E. Shevelev
LLC Scantronic Systems, Vvedenskogo kanala 7, 312, St. Petersburg 190013, Russia. E-mail: a [email protected], b [email protected].
This article describes the mathematical apparatus of original software that was developed by the authors for calculations of the radiation shielding and the fields of bremsstrahlung and scattered radiation of linear electron accelerators, which is adapted to a wide range of applications. The proposed program is based on the Monte Carlo method; this allows one to carry out a statistical calculation of the bremsstrahlung photon flux and to choose the trajectories of bremsstrahlung photons using a random number generator. The program is intended for the design of beam collimation systems, the calculation of complex radiation shielding profiles with labyrinths and tunnels, the optimization of the size of shielding barriers, assessment of the expected detection system response, and the prediction of the dose distribution in a cross section of a protected object. For verification of the reliability of calculations that are carried out by using this original software, the results of numerical experiments with the proposed code were compared to the results that were obtained using the Monte Carlo code MCNP version 4C3 for backscattered photons and forward transmitted bremsstrahlung photons with a boundary energy of 24 MeV for a number of materials and various thicknesses.
Keywords: bremsstrahlung, bremsstrahlung field, radiation shielding, linear electron accelerator, Monte Carlo
method, X-ray inspection systems, scattering, trajectory, photon, dose, transmission, albedo.
PACS: 28.41.Qb.
Received 17 November 2014.
English version: Moscow University Physics Bulletin 2(2015).
Сведения об авторах
1. Петрунин Владимир Иванович — докт. техн. наук; e-mail: [email protected].
2. Огородников Сергей Анатольевич — тел.: (812) 326-79-45, e-mail: [email protected].
3. Арлычев Михаил Анатольевич — e-mail: [email protected].
4. Шевелев Игорь Евгеньевич — e-mail: [email protected].
