Алгоритм расчета пространственного и спектрального распределений электронов, порожденных при взаимодействии проникающего излучения с объектами

( The Algorithm for Calculating Spatial and Spectral Electron Distributions Produced by Interaction of Penetrating Radiation with Objects
Preprint, Inst. Appl. Math., the Russian Academy of Science)

Введение.

Рассматривается задача о формировании электронных потоков внутри объектов, находящихся под воздействием гамма-излучения. Необходимость анализа пространственного распределения потоков релятивистских электронов, инжектируемых с внутренних поверхностей трехмерных объектов, возникает, например, при разработке расчетного аппарата для математического сопровождения экспериментов по облучению объектов ионизирующим излучением. В этих экспериментах измеряемой величиной является электромагнитное поле, генерируемое внутри облучаемого объекта вследствие движения электронов, порожденных при взаимодействии проникающего излучения с материалами оболочки объекта.

В рассматриваемых задачах основными типами взаимодействия фотонов излучения с конструкционными материалами объектов являются фотопоглощение, когерентное (Рэлеевское) и некогерентное (Комптоновское) рассеяния, рождение электрон-позитронных пар [1]. Все указанные процессы, кроме когерентного рассеяния, приводят к рождению высокоэнергетичных электронов, которые могут инжектироваться с внутренних поверхностей объекта, что приводит к возбуждению электромагнитного поля внутри. Для расчета этого поля необходимо корректное задание пространственного и энергетического распределений плотности тока инжектируемых электронов.

В настоящей работе предложен алгоритм для расчета плотности потока вылетающих внутрь объекта электронов для практически важных случаев, когда преобладающим механизмом взаимодействия проникающего излучения с веществом оболочки объекта является комптоновское рассеяние. При исследовании радиационного возбуждения электромагнитных полей в технологических объектах часто имеет место ситуация, когда толщина внешней оболочки такого объекта значительно меньше средней длины пробега фотона в соответствующем веществе. Построенный алгоритм в явном виде учитывает это соотношение, что позволило существенно повысить эффективность расчетов за счет пренебрежения маловероятными эффектами.

В основе разработанного алгоритма лежит использование модифицированного метода Монте-Карло как для расчета трансформации проникающего излучения в конструкционных материалах объектов [2], так и для моделирования электронных траекторий. Сконструированный алгоритм реализован в виде параллельного программного комплекса с использованием интерфейса MPI. В работе приведены примеры результатов расчетов пространственного и спектрального распределения электронов, полученных с помощью указанного программного комплекса.

Физическая модель.

Для описания процессов взаимодействия гамма-излучения с материалом объекта используется общепринятая физическая модель [1,3]. Предполагается, что основным типом взаимодействия фотонов с веществом, приводящим к образованию электронов, является комптоновское рассеяние. Рассматриваются практически важные ситуации, в которых энергии гамма-квантов и конструкционное вещество объекта таковы, что характерная толщина материала оболочки объекта много меньше длины пробега фотонов. В этой ситуации взаимодействие фотона с веществом оболочки является маловероятным событием, что позволяет ограничиться рассмотрением однократного рассеяния кванта внутри оболочки объекта, пренебрегая электронами, рожденными в результате вторичного и последующих столкновений.

Распределение длины траектории фотона в веществе определяется законом ослабления: , где  - коэффициент ослабления гамма-излучения в материале объекта. Вероятность взаимодействия фотона с тонким слоем вещества (,  - характерная толщина объекта) равна . Вероятность второго столкновения будет по порядку величины пропорциональна . Поэтому в рассматриваемых ситуациях учитывается только однократное взаимодействие фотонов внутри оболочки объекта.

В случае, когда взаимодействие фотона с веществом произошло, вероятность некогерентного рассеяния вычисляется по формуле , где  - коэффициент комптоновского рассеяния. Косинус полярного угла рассеяния фотона (рис.1)  – случайная величина, имеющая распределение Кляйна-Нишины [4]:

.

Используя закон сохранения энергии и импульса нетрудно получить формулы для «электронных» величин:

-   кинетическая энергия - , где  и  - энергия фотона в единицах  до и после взаимодействия соответственно;

-   полярный угол рождения электрона в системе координат, связанной с первоначальным движением фотона - ,  - косинус угла рассеяния фотона,  - электронный импульс в единицах ;

-   распределение длины пути электрона в веществе определяется законом ослабления: ,  - коэффициент ослабления электронного потока в материале объекта.

Подробное описание учитываемых физических процессов приведено, например, в [5].

Расчетный алгоритм.

Для моделирования процесса рождения и регистрации электронов (предполагается, что детекторная плоскость есть ) используется модифицированный метод Монте-Карло, в котором явно учитывается малый параметр задачи – вероятность взаимодействия гамма-квантов с тонким слоем вещества, что позволяет ограничиться учетом однократных столкновений.

При этом комптоновское рассеяние фотона рассматривается не как случайное событие, а как детерминированный процесс, причем рождающемуся комптоновскому электрону приписывается начальный «вес», равный вероятности указанного события. Это позволяет значительно уменьшить дисперсию результатов, а также упростить логику алгоритма.

Рассмотрим разработанный алгоритм расчета электронных потоков подробнее.

Пусть характерная толщина оболочки объекта равна . Длина отрезка , лежащего внутри материала, в направлении движения фотона -  (рис.1). Начальный «вес» электрона  вычисляется как произведение вероятности взаимодействия фотона внутри слоя вещества оболочки объекта толщиной   и вероятности того, что это взаимодействие будет комптоновским рассеянием , .

После вычисления начального веса электрона разыгрывается точка взаимодействия фотона с веществом оболочки. Для этого вначале определяется функция распределения  фотонов, рассеивающихся внутри оболочки,  - расстояние от точки входа кванта в объект до точки взаимодействия. Указанная функция имеет вид .

Затем, используя уравнение , где  - равномерно распределенная на  случайная величина [6], получаем формулу для розыгрыша координат точки взаимодействия фотона:

,

где  - направляющий вектор начального движения фотона,  – радиус-вектор точки входа фотона в материал объекта.

После определения координат точки взаимодействия кванта  осуществляется розыгрыш угла рассеяния  фотона в локальной системе координат, начало которой совпадает с найденной точкой , а ось  направлена вдоль первоначального движения фотона (рис.2).

Для определения указанного угла используется следующий алгоритм. Имея нормированное распределение Кляйна-Нишины , и равномерно распределенную в интервале (0,1) случайную величину ,  можно получить из уравнения (см., например, [6]):

Основную трудность при определении угла  представляет решение уравнения . Вычисление  для каждой фотонной траектории, учитывая, что число таких траекторий может достигать 10¹⁰ и более, требует огромных вычислительных затрат. Поэтому в алгоритме предусмотрено использование заранее подготовленных двумерных таблиц , построенных с заданной точностью на равномерных сетках по  и . При этом внимание уделяется как точности вычисления соответствующих интегралов, так и требованию малости изменения значений  в соседних ячейках таблиц. Методика построения указанных таблиц основана на численном интегрировании сеточных функций, построенных путем интерполяции экспериментальных данных [7], а также с помощью использования метода Ньютона определения значений функции, заданной неявно.

Применение этих таблиц позволяет свести вычисление косинуса полярного угла рассеяния фотона к нескольким арифметическим операциям определения номера элемента двумерного массива .

Замечание. Нетрудно построить аналогичные таблицы для непосредственного розыгрыша направления движения электрона без определения угла рассеяния фотона. Однако, использование «электронных» таблиц затруднительно в силу того, что угловое распределение рождающихся электронов обращается в ноль при  (т.е. при малых углах комптоновского рассеяния фотонов), что приводит к неограниченному возрастанию производной функции  и, как следствие, к заметной потере точности табличного определения электронного полярного угла в окрестности .

Отметим далее, что азимутальный угол  направления движения рассеянного фотона и угол  направления движения электрона связаны соотношением , поэтому для определения электронной траектории достаточно разыграть только один из этих углов. В разработанном алгоритме предусмотрен розыгрыш , поскольку фотонная траектория после рассеяния не рассматривается. Для однозначного определения начала отсчета угла  ось  выбирается таким образом, чтобы соответствующий орт лежал в плоскости , образованной осью  исходной системы координат и вектором, определяющим начальное движение фотона. Направление орта , выбирается так, чтобы при повороте плоскости  вокруг оси  на угол  направление орта оси  совпало с направлением орта оси .

Таким образом, для каждой фотонной истории разыгрываются три независимые случайные величины: расстояние  от точки входа кванта в объект до точки взаимодействия фотона с веществом оболочки объекта, косинус угла комптоновского рассеяния кванта  и азимутальный угол направления электронной траектории .

После определения полярного угла траектории рассеянного кванта и азимутального угла траектории электрона вычисляется энергия рассеянного фотона , косинус угла рождения  и кинетическая энергия электрона  (см. выше).

Направляющий вектор движения электрона в исходной системе координат , , матрица  перехода из локальной системы координат в исходную определяется двумя поворотами: сначала плоскость  поворачивается вокруг оси  на угол , а затем вокруг оси  на угол . В этом случае  имеет вид:

.

На следующем шаге алгоритма определяется точка пересечения траектории электрона с детекторной поверхностью: , , расстояние  от точки взаимодействия фотона до детекторной плоскости и вероятность достижения этой плоскости электроном .

Наконец, определяется «фотонный» вес регистрируемого электрона в соответствии с угловым  и спектральным  распределениями проникающего излучения источника, когда рассматриваемый фотон испускается источником в направлении  с энергией  (этот шаг алгоритма можно, очевидно, выполнять и в начале рассмотрения каждой фотонной траектории).

Для этого вводится разбиение исходного фотонного энергетического интервала: , и «спектральный» фотонный вес  вычисляется следующим образом . Аналогично вычисляется «пространственный» фотонный вес электрона . В случае, если используется плоский источник, для определения  вводится разбиение общей площади , с которой испускаются фотоны, на элементарные площадки  и .

Конечный вес электрона , пересекшего детекторную поверхность, определяется следующим образом: .

Плотность электронного потока  в расчете на один фотон в фиксированном детекторном пикселе вычисляется по формуле: ,  - площадь регистрирующего пикселя на детекторной плоскости,  - общее количество фотонных историй,  - количество пересечений электронных траекторий с рассматриваемым пикселем.

Часто требуется определить не поток, а число электронов , пересекающих данный элемент детекторной поверхности, приходящееся на один гамма-квант. Формула для расчета этой величины имеет вид: , причем  зависит от пространственного распределения излучаемых фотонов. В частности, для конусного источника  - это число фотонов, излученных источником в телесном угле, под которым «виден» рассматриваемый пиксель детекторной поверхности. В случае плоского источника, когда направление исходного излучения совпадает с нормалью к детекторной плоскости, , где  - площадь поверхности источника.

Тестирование алгоритма.

Для качественной оценки результатов работы сконструированного алгоритма использовались очевидные свойства симметрии, которым должно удовлетворять получаемое пространственное распределение электронных потоков. Так, например, для плоского источника фотонов продольная компонента плотности потока должна быть равномерно распределена на детекторной плоскости при облучении пластины вдоль ее нормали. При использовании конусного источника излучения -компонента потока должна совпадать на срезе вдоль оси  с -компонентой на срезе вдоль оси  и т.п.

Рассмотрим конфигурацию эксперимента, когда коллимированный источник гамма-излучения падает на пластину толщиной  вдоль ее нормали (рис.3). В этой ситуации удается построить решение задачи о формировании потока электронов, порожденных в результате однократного рассеяния фотонов, в квадратурах.

Плотность электронного потока   в расчете на один фотон будет в этом случае определяться формулой:

.

Здесь  - фотонный спектр;  - линейный коэффициент комптоновского рассеяния ();  - линейный коэффициент ослабления гамма-излучения ();  - нормированное угловое распределение комптоновских электронов, порожденных фотоном с энергией ; ;  - линейный коэффициент ослабления электронного потока ().

Угловое распределение комптоновских электронов определяется следующим образом. Вероятность рождения электрона (в случае, когда комптоновское рассеяние произошло), имеющего направление движения в диапазоне углов от 0 до  равна

, .

Эта вероятность должна, очевидно, совпадать с вероятностью рассеяния фотона в диапазоне углов от  до . Обозначим , причем функция  определяется из соотношений: , , , . Следовательно,

.

Отсюда .

Результаты сравнительных расчетов изображены на рис.4. На этом рисунке приведены графики зависимости плотности электронного потока (в относительных единицах) от расстояния соответствующего элемента детекторной поверхности от начала координат вдоль оси  в см. Сплошной линией показаны результаты расчетов с помощью разработанного алгоритма, кружками отмечены результаты тестовых расчетов. Расчеты проведены для монохроматического источника с энергией квантов равной 1 мэВ и алюминиевой пластины толщиной 1 см.

Анализ полученных результатов показал, что тестовые значения плотности электронного потока и соответствующие величины, полученные с использованием описанного алгоритма, совпадают в рамках статистической погрешности разработанной методики.

Результаты расчетов для плоского и конусного источников.

Для анализа эффективности разработанного алгоритма и оценки качества результатов были проведены расчеты электронных потоков, образующихся при облучении пластины проникающим излучением. В качестве примеров приведем результаты расчетов для квадратной железной пластины толщиной 2мм и двух источников тормозного излучения, имеющего типичный энергетический спектр. В расчетах использовались точечный источник, имеющий равномерное по телесному углу пространственное распределение излучения, и плоский источник площадью  с равномерным по площади распределением излучения.

Спектр использованного тормозного излучения изображен на рис.5.

Рассмотрим вначале результаты расчетов с использованием конусного источника проникающего излучения. Геометрия вычислительного эксперимента изображена на рис.6.

Расчеты пространственного распределения компонент электронного потока проводились на многопроцессорной вычислительной системе МВС-1000М. Для достижения уровня статистической погрешности результатов менее 0.5% учитывалось  фотонных траекторий. Время расчета одного варианта составило около часа при использовании 200 процессоров.

На рис.7 приведены графики полученных пространственных распределений декартовых компонент плотности потока комптоновских электронов в относительных единицах. На этой иллюстрации сплошной линией изображена зависимость -компоненты плотности потока электронов от номера детектора вдоль оси . Точками отмечена зависимость -компоненты потока вдоль оси . Третья компонента плотности потока практически равна нулю вдоль указанного направления. Пространственные распределения вертикальной () компоненты  и горизонтальной  являются аксиально-симметричными. В точке пересечения оси симметрии конуса и детекторной плоскости  - максимальна, в то время как  обращается в ноль. Высота расположения источника подобрана таким образом, чтобы угол между осью  и образующей конуса, проходящей через абсциссу края пластины, был равен . Поэтому вертикальная и горизонтальная компоненты плотности электронного потока совпадают по абсолютной величине в токах детекторной плоскости, соответствующих точкам пересечения осей  и  с краями пластины.

Указанные выше свойства симметрии полученных результатов соответствуют, очевидно, априорным оценкам свойств симметрии искомого пространственного распределения компонент плотности электронного потока.

Вторая серия расчетов была проведена для плоского источника, расположенного параллельно детекторной плоскости. В проведенных расчетах использовалось проникающее излучение с тем же спектром (см. рис.5) и такая же, что и в случае конусного источника, железная пластина.

Отметим, во-первых, что полученные результаты удовлетворяют с хорошей точностью очевидным свойствам симметрии задачи. В частности, горизонтальная компонента  плотности электронного потока близка к нулю практически на всей детекторной плоскости, за исключением малой окрестности границ теневой области, образованной проекцией площади источника на детекторную поверхность. Кроме того, вертикальная составляющая  имеет постоянное значение в пределах указанной теневой области для любого энергетического диапазона рассчитанного электронного спектра.

На рис.8 приведен групповой спектр  вертикальной компоненты полученного электронного потока. По оси абсцисс отложена энергия в мэВ, по оси ординат – групповой спектр электронов в относительных единицах. В силу того, что средняя длина свободного пробега электрона для большей части комптоновских электронов превышает толщину пластины, электронный спектр близок по форме к гамма-спектру.

На рис.9 изображена -компонента плотности рассчитанного электронного потока, проинтегрированная по спектру. Пространственное распределение указанной компоненты изображено для наглядности в виде поверхности над детекторной плоскостью, причем по осям  и

 отложены соответствующие номера элементов детекторной матрицы, а по оси  - значение  в относительных единицах.

Приведенные в качестве примера иллюстрации позволяют сделать вывод о передаче основных априорных свойств симметрии задачи, а также о достаточно высокой эффективности предложенного алгоритма, позволяющего рассчитывать с достаточной точностью электронные потоки, порожденные при взаимодействии проникающего излучения с веществом. При этом соответствующие расчеты требуют приемлемых вычислительных затрат на современной многопроцессорной технике.

Заключение.

Разработан алгоритм расчета пространственного распределения потоков релятивистских электронов, инжектируемых с внутренних поверхностей трехмерных объектов в экспериментах по облучению этих объектов -излучением. Указанный алгоритм реализован в виде параллельного программного модуля для проведения расчетов на многопроцессорной вычислительной системе МВС-1000М.

Предложенный алгоритм в явном виде учитывает малый параметр задачи – вероятность комптоновского рассеяния -квантов выбранного энергетического диапазона внутри оболочки объекта, что позволяет значительно увеличить число используемых в расчете фотонных историй для уменьшения статистической погрешности результатов без увеличения вычислительных затрат. Кроме того, указанный подход дает возможность значительно упростить логическую часть алгоритма, что повышает надежность и быстродействие программного модуля. Построен трехмерный аналитический тест для проверки точности результатов расчетов. Проведено сравнение полученных результатов с указанным тестом, которое показало их совпадение в рамках статистической погрешности.

Литература.

1.            У.Фано, Л.Спенсер, М.Бергер. «Перенос гамма-излучения», Госатомиздат, Москва, 1963.

2.            M.Zukhovsky, S.Podoliako, G.-R.Tillack, and C.Bellon: Monte Carlo Simulation of Photon Transport Coupled to CAD Object Description. In: D. O. Thompson and D. Chimenti (Eds.): Review of Progress in Quantitative Nondestructive Evaluation, Vol. 23, AIP, Melville, 2004, pp. 515-521

3.            Н.Г.Гусев, Л.Р.Кимель, В.П.Машкович, Б.Г.Пологих, А.П.Суворов. Защита от ионизирующих излучений, Атомиздат, Москва, 1969.

4.            А.И.Ахиезер, В.Б.Берестецкий Квантовая электродинамика, Издательство «Наука», Москва, 1969.

5.            С.В.Подоляко, Е.Г.Лукьянова. Численное моделирование трансформации рентгеновского излучения в объектах с учетом влияния форм-факторов на угловое распределение фотонов. Препринт Института прикладной математики им. М.В.Келдыша РАН, №6, 2004.

6.            И.М.Соболь. Метод Монте-Карло. «Наука», Главная редакция физико-математической литературы, Москва, 1972.

7.            J.F.Briesmeister (ed.): MCNP–A General Monte Carlo N–Particle Transport Code. LANL Report LA–13709–M, Los Alamos, 2000.