автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Прецизионные методы Монте-Карло для расчета транспорта электронов

На правах рукописи УДК 519.85+614.876+621.039.5

БЕЛОУСОВ ВИКТОР ИВАНОВИЧ

ПРЕЦИЗИОННЫЕ МЕТОДЫ МОНТЕ-КАРЛО ДЛЯ РАСЧЕТА ТРАНСПОРТА ЭЛЕКТРОНОВ

13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

1 2 АПР 2012

Обнинск-2012

005018436

Работа выполнена в Обнинском институте атомной энергетики - филиале Национального Исследовательского Ядерного Университета «МИФИ».

Научный руководитель:

Официальные оппоненты:

Ведущая организация:

доктор физико-математических наук, профессор ИАТЭ НИЯУ МИФИ Кураченко Юрий Александрович

доктор физико-математических наук, главный научный сотрудник ГНЦРФФЭИ.

Коробейников Валерий Васильевич

доктор физико-математических наук, главный научный сотрудник РНЦ «Курчатовский институт». Гуревич Михаил Исаевич

Научно-исследовательский институт ядерной физики имени Д. В. Скобельцына Московского государственного университета имени М. В. Ломоносова.

' Защита диссертации состоится « 27 » апреля 2012 г. в 10 часов на заседании диссертационного совета Д 201.003.01 в ГНЦ РФ- Физико-энергетическом институте имени А. И. Лейпунекого по адресу: 249033, Калужская область, г. Обнинск, пл. Бондаренко, д. 1.

С диссертацией можно ознакомиться в библиотеке ГНЦ РФ-ФЭИ.

Автореферат разослан «Н» марта 2012 г.

Ученый секретарь диссертационного совета, доктор технических наук

^^Г Т.н. Верещагина

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ.

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

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

Цель работы состояла в дальнейшем развитии комплекса программ BRAND, как прецизионного инструмента численного решения задач переноса электронов. В рамках данного ПК, используя полученный модуль подготовки константной информации по файлам оцененных ядерных данных библиотек EEDL [1] и EADL [2] формата ENDF-6 [3], были разработаны и программно реализованы алгоритмы моделирования процесса переноса электронов. Метод моделирования индивидуальных соударений исключает дополнительное осреднение потока частиц [4-6], что позволяет применить классическую модель метода Монте-Карло. Неопределенность в оценке получаемого численного решения сводится к неопределенности константных данных, т.к. фактически «напрямую» используется информация из файлов библиотек оцененных ядерных данных. Константный модуль ПК BRAND является уникальным, т.к. в нем удалось объединить всю необходимую информацию о взаимодействиях электронов с атомами вещества при условии бинарных взаимодействий. Проверка полученного константного модуля и алгоритмов моделирования взаимодействий электрона с материалом были сделаны благодаря применению полуаналитического метода (ПА-метод) [7-9], который дает возможность имитировать энергопотери частицы в веществе без моделирования длины свободного пробега, что существенно уменьшает погрешность вычислений, и фактически решение задачи переноса представляется в аналитическом виде. Научная новизна настоящей работы заключается в следующем:

1. Разработана трехмерная математическая модель имитации задач переноса электронов -метод индивидуальных соударений (.МИС), который является разновидностью классического моделирования методом Монте-Карло.

2. Созданы и программно реализованы алгоритмы моделирования процессов рассеяния электронов, использующие информацию файлов оцененных ядерных данных библиотек EEDL и EADL формата ENDF-6.

3. Разработана новая версия ПК BRAND для решения задач переноса электронов с применением МИС и ПА-метода.

Практическая значимость заключается в применении предложенных теоретических методов для расчета транспорта электронов в рамках ПК BRAND, при этом созданный

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

Личный вклад автора:

1. Построены алгоритмы моделирования переноса электронов ПК BRAND с применением прецизионных методов Монте-Карло: МИС и ПА-метод.

2. Разработаны алгоритмы моделирования рассеяния электронов с использованием константной информации непосредственно из файлов оцененных данных формата ENDF-6.

3. Полученные алгоритмы реализованы и программно интегрированы в ПК BRAND.

4. Проведены расчеты, численно имитирующие транспорт электронов при условии бинарных взаимодействий в различных средах.

Результаты работы, выносимые на защиту

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

Приводятся примеры решения задач переноса электронов, в том числе задач глубокого проникновения частиц и задача внутрисосудистой радиотерапии [10-12].

Основные положения диссертации, выносимые на защиту:

1. Математическая модель - моделирование индивидуальных соударений электронов с атомами вещества или МИС.

2. Численные методы - прецизионные методы Монте-Карло, модифицированные для решения задач переноса электронов -МИС и ПА-метод.

3. Комплекс программ - ПК BRAND, приспособленный для решения задач транспорта электронов.

4. Результаты решения задач транспорта электронов с помощью ПК BRAND и сравнение с аналогичными ПК.

Апробация работы

Основные результаты работы опубликованы в виде 10 статей (4 статьи в журналах из списка ВАК). По материалам исследований были представлены 11 докладов на российских и международных научных конференциях. Список публикаций и докладов приведен ниже.

Структура и объем диссертации

Диссертация состоит из введения, четырех глав, заключения, списка литературы из 59 наименований, списка таблиц и списка рисунков. Общий объем работы составляет 178 страниц, включая 43 рисунка, 4 таблицы и 9 приложений.

СОДЕРЖАНИЕ РАБОТЫ.

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

В первой главе рассматриваются основные теоретические вопросы задач теории переноса излучений. Описываются базовые принципы решения уравнения переноса методом Монте-Карло. Дано теоретическое описание применения метода Монте-Карло в задачах переноса излучения, в частности задачи переноса электронов. Отмечается преимущество реализации метода Монте-Карло для решения уравнения переноса электронов в рамках дозиметрических задач. Приведено детальное описание основных методов решеиия задачи переноса электронов, используемых на сегодняшний момент в программных комплексах MCNP [4] и PENELOPE [5]. Дано алгоритмическое описание методов. Определены недостатки алгоритмов решения задач переноса электронов. Сформулирована концепция разновидности метода Монте-Карло - метода моделирования индивидуальных соударений частиц, специально разработанного для обработки информации о взаимодействии электронов с атомами вещества, которая содержится в файлах библиотек оцененных ядерных данных. Предлагается сравнение предложенного метода с аналогами. Обосновывается, почему при моделировании задач переноса электронов предлагаемый вероятностный метод дает более точное решение.

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

Рассматриваемые дозиметрические задачи формируют сложную трехмерную геометрию системы. Более того, конфигурация решаемой проблемы существенно меняется в зависимости от заданных нуклидов и их свойств. Поэтому использование метода Монте-Карло является единственным приемлемым решением проблемы, т.к. указанный метод не только одинаково эффективно работает в условиях геометрий произвольной сложности, но также приспособлен к детальному учету всей имеющейся информации, т.е. обеспечивает наилучшую точность, обусловленную точностью исходных данных. Именно эта идеология применения «истинного» метода Монте-Карло при решении уравнения переноса излучения и реализована в ПК BRAND.

Краткими итогами первой главы являются следующие основные положения.

Дано теоретическое описание применения метода Монте-Карло в задачах переноса излучения, в частности задачи переноса электронов.

Приведено детальное описание основных методов решения задачи переноса электронов. Дано алгоритмическое описание методов. Определены недостатки методов и алгоритмов решения задач переноса электронов.

Сформулирована концепция разновидности метода Монте-Карло - МИС.

Во второй главе предлагается к рассмотрению подробное описание алгоритмов моделирования взаимодействия электронов с атомами вещества. Информация о взаимодействии электрона с заданным материалом используется напрямую из файлов оцененных ядерных данных формата ENDF-6 библиотек EEDL и EADL.

При выполнении исследования использовались текстовые файлы оцененных данных формата ENDF-6, в частности файл вторичных распределений оцененных фотонно-электронных ядерных данных (Файл 26) и файл данных атомной релаксации (Файл 28). Файл 26 содержит данные для описания следующих реакций электрона: упругое взаимодействие, тормозное излучение, возбуждение атома и электроионизация атома. Поэтому были разработаны четыре основных алгоритма моделирования характеристик реакций электрона согласно формализмам формата ENDF-6, а также алгоритм моделирования хода реакции электроионизации с последующей релаксацией атома

Описания получения некоторых сечений можно взять из [13], раздел «Взаимодействие заряженных частиц с веществом». На основе теоретических выкладок раздела «Electron Interactions», взятых из [4] и описаний Файлов 26 и 28, была проведена разработка алгоритмов и программ по моделированию взаимодействий электрона. Краткое описание алгоритмов представлено далее.

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

Для моделирования косинуса угла отклонения вылетающего электрона известно несколько плотностей углового распределения^, р), зависящих от энергии налетающего электрона и дающих, при определенном Ео, некоторую функцию плотности углового распределения, зависящую от косинуса угла вылетающего электрона (ц = cos(fl)). Так как функции плотностей распределения заданы не непрерывно, а дискретным набором энергий, то в начале определяются энергии, для которых заданы две функции распределения плотностей. Между ними находится заданная Е0 т.е. Et <Ео<Ег. Получаем две функции/5(£ь р) и Д(Е2, д). Известно, что между ними лежит искомая функция/^ ц), т.е.Д(Еи р) И) <МД1у\1).

После того как функции плотности углового распределения определены, необходимо^ найти Цо в зависимости от смоделированного случайного числа F, причем F = у ~ [/(0,1). Основным свойством данной случайной величины является следующее

двух случаев, т.е. для двух функций плотности распределения^,, р) и Д(Е2, ц) и находим два решения ц, и ц2 соответственно. Известно, что ц, < р0 < р2. Также известно, уо точки (£,, р,), (Е2, р2) и (Е0, ро) лежат на одной прямой. Поэтому находим ц0 по формуле

соотношение

Решаем эту задачу для

(1)

Соотношение (1) следует из решения системы уравнений:

ц0=а-Е0 + Ь Гц0-ц, =а-(£:0-К1)

В результате получаем, что

а отсюда искомое выражение (1).

Функция плотности углового распределения представлена дискретным образом на отрезке поэтому оставшийся до единицы отрезок [ц',1] заполняется аналитиче-

ской функцией плотности:

Л(£.м) = Т—(4) (л+1-ц)

Параметр Л можно выразить явным образом через тр

Л = /(т1 + 1-ц')\ (5)

где /=/4(Я,ц') - значение функции плотности углового распределения в точке д'. В свою очередь параметр т) [4] вычисляется выражением:

1 ( а-т-с | _у 4 = r-hrm—\-ZH-

2 I 0.885 -р

где

I.I3+3.76-(a-Z/ß)2-^-lj.

(6)

(7)

а а- константа тонкой структуры, т - масса покоя электрона, р - импульс электрона, v - скорость электрона, а Z- атомный номер.

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

В и. 2.4.1 изложено моделирование энергии порожденного фотона, для которого аналогично упругому взаимодействию заданы несколько плотностей углового распределения £р), зависящих от энергии налетающего электрона и дающих, при определенном Ео, некоторую функцию плотности распределения энергии, зависящую только от энергии вылетающего фотона (£р). Аналогично упругому взаимодействию переписываем соотношения, заменяя ц на Ер. Необходимо учитывать, что £,<Еа< Ег Получаем две функции /,(£ь Е^ и Д(Е2, Ер), между которых лежит искомая функция Д(Ео, Ef\т.е./5(£,, £р) £„) </?(£2, £,).

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

F = y~ U(0, 1). Причем F= ] fK(Eü,Ep)dEp = \ f№,Ep)dEp=\ ф1ъЕр)с1Е, . Нахо-

Ooo дим решение этого уравнения для двух случаев, т.е. для двух функций плотности распределения Е^)

и №. Ер) я получаем значения двух параметров E^i и £'р2 соответственно. Известно, что £р, < Еро < Ер2. Известно, что точки (Ei.-Epi), (Е2,Ер2) и (Ео, £ро) лежат на одной прямой. Поэтому находим Е& по формуле

--Е _ „-, которая следует из решения системы уравнений:

Е„ = а-Е, + Ь Егг = а-Е2+Ь

к = „•(£„-£,)

Заметим, что косинус направления движения фотона цр определяется равенством:

где у~(/(0,1),а/? берется из соотношений (7).

В п. 2.4.2 дано описание моделирования энергии вылетающего электрона, для

которого заданы энергии, вычитаемые из исходной энергии налетающего электрона (Е<>). Для определения вычитаемой энергии производится поиск двух соседних энергий с Еа, т.е. и Е2, по заданной дискретной функции.

Вычисляем искомое /о =.Д£о) по следующей формуле аналогичной (3)

= / + -(у;-/). Рассуждения здесь аналогичны рассуждениям, для поиска ц0-В итоге Ее~Ео -/о-

В разделе 2.5 изложена методика моделирования реакции возбуждения атома. При этом типе взаимодействия электрон теряет часть своей энергии и передает ее остаточному атому. Направление движения электрона при этом не меняется.

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

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

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

Алгоритм моделирования энергии «освобожденного» электрона (Ее) повторяет п. 2.4.1 предыдущего раздела, для порожденного фотона, поэтому приведем основные соотношения для моделирования:

Е^Е^Е, ЕЛ<.Еа4Е.г

■•г

(10)

(Е0,Б,)<1Е.=1А = ] А

р=

Аналогично (1) находим искомое Е,„ =

Энергия рассеянного электрона определяется по формуле:

(И)

Еш. = Е„ - Еш - Е,0,

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

В разделе 2.6 изложены общие принципы моделирования релаксации атома. Атом может переходить в состояние релаксации вследствие нескольких видов фотонно-электронного взаимодействия. При этом происходит передача некоторого количества энергии атому, которое в свою очередь распределяется по электронным уровням атома, выделяя радиационные фотоны или «выбиваемые» с более высоких уровней электроны (Файл 28).

Итак, подведем краткие итоги второй главы.

Глава содержит описание алгоритмов извлечения и применения информации файлов оцененных ядерных данных.

Предварительное описание файлов оцененных ядерных данных с последующим извлечением исходной информации из этих файлов позволило:

• дать подробное описание алгоритмов для расчета параметров движения частицы во всех реакциях взаимодействия электрона с атомом вещества;

• построить алгоритм протекания релаксации атома с последующим порождением дополнительных частиц моделирования переноса излучения.

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

Третья глава содержит анализ применения ПА-метода при расчете энерго-угловых спектров прямого и отраженного излучения электронов для однородного бесконечного барьера заданной толщины с1 (допустим случай й - со).

Применение в задачах переноса излучения аналогового метода Монте-Карло обычно требует значительного объема вычислительной работы, т.к. дисперсия метода велика. Одним из возможных путей, приводящих к уменьшению величины (/ - время реализации одной истории на компьютере, дисперсия на одну историю), является использование ПА-метода Монте-Карло, предложенного в [7, 8].

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

В разделе 3.6 кратко изложена методика получения рекуррентных соотношений для барьера заданной толщины. Пусть рк(г) - вероятность того, что частица находится в точке г перед ¿-ым взаимодействием. Тогда вероятность провзаимодействовать в элементе (к около точки г есть СкРк^)еЬ, где Ск X, - полное сечение к-го

взаимодействия, ч>к - косинус угла между направлением к-го рассеяния О, и осью 02, нормальной к поверхности барьера. Следуя [7], получаем

(12)

где Sk - сечение рассеяния ¿-го взаимодействия, т.е. Як/Ек есть вероятность не поглотиться на к-ом взаимодействии. Интегрирование ведется на отрезке [0, г] для случая > 0 или на отрезке [г,^] для случая < 0. Нетрудно показать, что

А(«) = 1 + (13)

»¡>0 »¡<0

причем = 1. Суммирование ведется по соответствующим / от 1 до к. Интегрируя (12) с учетом (13), получаем следующие рекуррентные соотношения для случая и^, > 0

■УЛ..

4»; =

ЧК^-с,)'

¡ = Ък

ТА,

Я*

(14)

и для случая и^, < 0

Л*и - ~

УА.

КК^-с,)*

/=и

=- 5Х,

„-с.-'

-2Х>

(15)

Итак, движение частицы организуется следующим образом: знерго-угловая часть плотности перехода моделируется МИС, по полученным значениям Й, и вычисляются необходимые величины для формул (14) или (15), затем находятся коэффициенты разложения (13) для р^г). Моделирование пробега частицы заменяется аналитическим интегрированием по всевозможным пробегам.

В качестве оценки энергетического спектра альбедо можно использовать оценку

(16)

о, *<0

гДе х(') = ||' ¡¿о ~ФункцияХевисайда,/-число рассеяний.

Оценкой энергетического спектра по заданному направлению Й„ будет

1 = Ёл(Й4->П0)-^(0), (17)

»=1

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

В разделе 3.7 дается краткий вывод рекуррентных соотношений дня полупространства. Если геометрия задачи - полупространство, то рекуррентные формулы существенно упрощаются. При интегрировании в (12) для случая <0, т.е. на полуинтервале [г,00), получается, что коэффициент при ехр(С4„-2) тождественно равен нулю. Таким образом, аналог (13) здесь есть

»1>0

Интегрируя (1) с учетом (7), получаем

УЛ.,

А+им = ~х ("ш) • 2 А*и

Соответственно изменятся и оценки

S = (20)

Отметим, что дисперсия tj не зависит от угла между осью OZ и направлением О,,, что позволяет оценивать угловой спектр отраженного излучения при wD близких к нулю. Это свойство оценки т| аналогично свойству независимости дисперсии % от толщины барьера d (см. [9]).

В разделе 3.8 разъясняются проблемы вычисления коэффициентов через рекуррентные соотношения, а также предлагаются варианты их решения. Значения коэффициентов Ct = It/w„ практически равны друг другу, вследствие очень больших значений сечений и малых изменений углов. Существует необходимость сохранения заданной точности при вычислении коэффициентов, используемых затем в рекуррентных соотношениях и при вычислении оценок. Надо заметить, что точность зависит от максимального количества взаимодействий электрона с атомами вещества.

Коэффициенты рекуррентных соотношений, вычисленные с помощью стандартных вещественных типов данных, не могут считаться истинными, так как многократно возрастающая с каждой итерацией погрешность выходит за пределы допустимой для стандартных вещественных типов. Чтобы получить значения коэффициентов, надлежит использовать библиотеку многократно увеличенной точности GMP (GNU MP) [14]. Библиотека является свободной в использовании и распространяется на бесплатной основе с соблюдением авторских прав копирования и распространения.

В п. 3.8.1 дается краткое описание библиотеки GMP. GNU Multiple Precision (Многократно увеличенная Точность) - это портативная библиотека, написанная на языке Си под арифметику с произвольной точностью для целых, рациональных и вещественных чисел. Цель создания библиотеки - обеспечить возможную наиболее быструю арифметику для всех приложений, которым необходима более высокая точность, чем та, что поддерживается стандартными типами Си.

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

Итак, процесс моделирования переноса электронов в условиях барьерной геометрии разбивается натри последовательных этапа:

1. моделирование энерго-угловых характеристик переноса электрона. На этом этапе применяется МИС\

2. моделирование характеристик траектории переноса электрона через рекуррентные соотношения. На этом этапе применяется ПА-метод;

3. вычисление оценок энергетического спектра.

Подведем краткие итоги третьей главы.

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

Проблемы вычисления коэффициентов потребовали использования возможностей библиотек сохранения «бесконечной» точности GMP и MPFR [15].

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

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

первичных и вторичных электронов, а также энергетического спектра электронов прошедших через барьер и отраженных от барьера, т.е. альбедо.

Моделирование процесса переноса излучения электронов в кремнии. Пример с простои геометрией был взят из эксперимента [10,11]: оценка потоков моноэнергетических электронов для кремния. Исходные данные эксперимента следующие: барьер из кремния толщиной 0,048085 см, точечный источник находится в начале координат оценивается спектр электронов с начальной энергией 2,43 МэВ после прохождения через пластину. Количество дискретных разбиений по энергии - 600. Результаты вычислений можно видеть на рис. 1. Видно хорошее совпадение результатов расчетов и эксперимента (статистическая погрешность расчета - менее 1 %).

15,0-

-„ 12,5 аз

I 10,0 я

а 7-5

с

►в 5.0.

2,5

данные эксперимента [11] - моделирование BRAND

0,0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 1,8 2,0 2,2 2,4 Энергия электрона, МэВ

Рис. 1. Оценка плотности потока электронов с начальной энергией 2,43 МэВ после прохождения через кремниевую пластину толщиной 0,48085 мм

Моделирование процесса переноса излучения элекггронов в алюминии. Следующий пример был взят из эксперимента [10,12]: оценка мощности дозы моноэнергетических электронов в алюминии. Исходные данные эксперимента следующие: барьер из алюминия, точечный источник находится в начале координат, оценивается распределение мощности дозы электронов с начальной энергией 1 МэВ при прохождении данной пластины. Результаты вычислений можно видеть на рис. 2. Также на рисунке проведено сопоставление с расчетами программы PENELOPE [5] по данным того же эксперимента. Моделирование решения данной задачи с использованием программы PENELOPE проведено самостоятельно.

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

I.

1.0

S 0.5

. 1 —,--р-—-- • экспернмент [ 12]

/ BR/ VND

/

/

ч

—— ■ .-1 —

0,175

0,200

0,000 0.025 0.050 0,075 0,100 0,125

Глубина, мм

Рис. 2. Оценка распределения мощности дозы по глубине проникновения в алюминиевую пластину электронов с начальной энергией 1 МэВ

1иПр°145Ра1,»,ВеНН0е Р«пр«диение поглощенной энергии для изотопов "*Re,

Но, Dy, *Y. Для получения распределения поглощенной дозы [16,17] рассматривалась вспомогательная задача - расчет потерь энергии в воде в зависимости от расстояния до точечного источника излучения. Моделирование решения данной задачи проводилось независимо с помощью двух программных комплексов - BRAND и MCNP. Полученные результаты сравнивались с соотношением Левингера [18], описывающим распределение поглощенной дозы для точечного изотропного ^-источника в неограниченной однородной тканеэквивалентной среде.

Для рассматриваемой задачи и выбранных изотопов выполняются дополнительные условия:

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

• внесение р-источника в среду не должно нарушать ее гомогенности, т.е. плотность и атомный номер поглотителя после внесения р-источника в среду должны изменяться незначительно;

• р-излучатель должен иметь простой эмиссионный спектр.

1«ПоЛ?б5еНН9^ РезУльтаты расчета потерь поглощенной энергии для изотопов IMRe,

Но, Dy, ^ представлены на рис. 3-6.

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

9 8 7. 6

5 4 3 2 1 О

0,00 0,02 0,04 0,06 0,08 0,10 0,12 0,14 0,16 0,18 0,20 Расстояние от источника, см

Рис. 3. Поглощенная доза от точечного источника "*Re Доза, 10'7, Гр/распад

7| 1 j ; | 1111-Г"!-

........___ - -функция Левингера [18]

6--—|--1--------расчёт MCNP -

.......1............|............ .................................. .расчёт BRAND

з-—j—--j------j--1--:---—

2----:------j--1------

о-------i--.----------

0,00 0,02 0,04 0,06 0,08 0,10 0,12 0,14 0,16 0,18 0,20 Расстояние от источника, см Рис. 4. Поглощенная доза от точечного источника ,44Но

Доза, 10'7, Гр/распад

... ........;...... ........j........ ! ..... ...... — функция Левингера [18] -расчёт MCNP --расчёт BRAND

...... 1 ......i......

........j......- ........j...

i ......;..... ........ - ........i..... ................ ...... ........ .....1—

......... ......... .........i- .... ....... ................ .......... ......;.........

....;...... .......t...... .....i........

v..... ....... .......i...... ....... ....... ....... .........!...... ......!.......

\i ........г ...... ........ !

........ \ ...... i

! :

Доза, 10', Гр/распад

-......;.......... ........ ........ ........ Г 1 >11 - функция Левингера [18] - расчёт MCNP — - »MinirBf RR ДХГП

......Л

......\ ........t...... ......i...... ...... ...... ...... '.........1........ ........1........ ......

.....\ ........!..... ....... .... .....j......... i ........;....... ......

..............* .........i.......... ....... ....... .......... ! ........f........ """Г....... .......

\ i i

........ ....... ...... i .......t........

10 9 8 7 6

5

4

3 2 1 О

0,00 0,02 0,04 0,06 0,08 0,10 0,12 0,14 0,16 0,18 0,20 Расстояние от источника, см Рис. 5. Поглощенная доза от точечного источника ttsDy

Доза, 10"7, Гр/распад

8Н-Щ--i--1--i--i-1-1-¡-1-:-1-:--

......~ ~ ! i -функция Лсвилгсра [18]

7--:-----i--i--i-----расчёт MCNP —

..........¡...........................................- ..............j......... .расчет BRAND

6--—I---;--:--!--i---1--;--j—-

5--Г--¡--;--i--j---1--г----

4--1—|--1--~--j--L---j--¡--j--

3--j---!--;--1--i---1--;--1--

0,00 0,02 0,04 0,06 0,08 0,10 0,12 0,14 0,16 0,18 0,20 Расстояние от источника, см Рис. 6. Поглощенная доза от точечного источника "Y

Доза, 10'7, Гр/распад

р -1-1-;-1-:-!• —- -функция Левингера [18] — - расчёт MCNP -пяг.црт RR AMD

\ .............

i .............. .... ¡...

' 1 i ......I . . ■ :' ■ j

......• 1 1...... .......;.....

i 1 1 ; ■■■!...... í

.....[......| ; ! ........... ........i........

..... j- \......... ......|..... j

Л ; .. к

....... ; .....i......

Доза, 10", Гр/распад

.....H= ......j........ ..........j........ i i '.■)■'.'■■ — функция Левингера [18]

Г — расчс -г BRA] VD

................... .......j......... .........

.....4 -' -■;...... .......i-....... .......i-...... ........;......

.......|........ ..........I. .. .........j........ .....|

\ ...... i _ ........f......

1 • V -!

\ H'' i ■■-! -

Задача внутрисосудистой радиотерапии. Следующий пример с реальной геометрией [19-23] относится к области задач внутрисосудистой терапии, связанных с радиохирургией. Практический интерес данной задачи сопряжен с использованием полученных результатов при исследовании параметров оборудования, предназначенного для лечения злокачественных опухолей и различных артериальных заболеваний, например, таких как артериальный стеноз (закупорка сосудов).

Исследуемая задача упрощенно формулируется следующим образом: катетер, который считается единым целым без механизма центрирования в сосуде, предположительно находится внутри артерии, заполненной кровью. Подробное описание геометрии исследуемой модели можно найти в [19], откуда можно узнать размеры каждого компонента, которые даны на соответствующих рисунках, в то время как композицию материалов и массовую плотность можно найти в прилагаемой таблице. Основные характеристики источника 32Р следующие: длина источника - 27 мм, активность - 7 ГБк, максимальная энергия частиц -1,71 МэВ.

Исследование проводилось с целью решения нескольких проблем поставленной задачи: доза, получаемая на глубине проникновения 1 мм и 0,5 мм от внутренних стенок артерии при отсутствии тромба на стенке артерии или при наличии тромбов различной плотности. Подробнее об исследуемых вопросах рассматриваемой задачи и полученных смоделированных решениях можно узнать из итогового отчета [24]. Результаты сравнительного моделирования представлены на рис. 7-10. При моделировании использовались следующие параметры: количество историй - 10, количество дискретных разбиений по расстоянию от центра источника- 100.

Доза, 10, Гр/час

-2,0 -1,5 -1,0 -0,5 0,0 0,5 1,0 1,5 2,0 Рис. 7. Распределение поглощенной дозы излучения электронов вдоль оси источника на расстоянии 1 мм от внутренней стенки артерии с тромбом и без тромба

Доза, 102, Гр / час

•В— данные [22] тромб-50%

расчет BRAND тромб-50% Ф данные [21] тромб-100% А- данные [22] тромб-100%

.........расчет BRAND тромб-100%

Э- данные [22] тромб-150% —- расчёт BRAND тромб-150%

5н»еее«е-»ее( >еввв« юеве* Расстояние от центра источника, см -1,0 -0,5 0,0

Рис. 8. Распределение поглощенной дозы излучения электронов вдоль оси источника на расстоянии I мм от внутренней стенки артерии с тромбом разной плотности

Лоза. 10*. Го/

час

-Г1—.-г

-2,0 -1,5 -1,0 -0,5 0,0 0,5 Рис. 9. Распределение поглощенной дозы излучения электронов вдоль оси источника на расстоянии 0,5 мм от внутренней стенки артерии с тромбом и без него

Доза, 102, Гр / час

—В—данные [22] тромб-50%

-расчет BRAND тромб-50%

данные [22] тромб-100%

-----------расчет BRAND тромб-100%

® данные [21] тромб-100% - в- данные [22] тромб-150% расчет BRAND тромб-150%

Ф "О.....Ф;.....QÏ

Ф "!®......о- е-

веек >&©iee-i »оооск >

Расстояние от цапра источника, см

_ -1,0 -0,5 0,0 о;5 1.0 1,5

1а- Распределение поглощенной дозы излучения электронов вдоль оси источника на расстоянии 0,5 мм от внутренней стенки артерии с тромбом разной плотности г

35 Доза, 10* Гр/час

3025-

расстояние: 1,0мм расстояние: 0,5ми

-без тромба---без тромба

.....тромб-100% -----тромб-100%

-........тромб-50% --------тромб-50%

■ тромб-150% -----тромб-150%

» „ D -1'5 -1-0 -0-5

Распределение поглощенной дозы излучения электронов вдоль оси источника на расстояниях 1,0 и 0,5 мм от внутренней стенки артерии с тромбом разной плотности

Моделируемое решение достаточно точно аппроксимирует авторское решение задачи [21, 22]. Погрешность полученных результатов объясняется статистической погрешностью моделирования и не превышает 2 % на интервале не более 1см от центра источника, заданного вдоль оси. Косвенным подтверждением достоверности полученного решения может служить геометрическая несимметричность данной модели, что видно из предложенных для сравнения рисунков. Моделирование BRAND сохраняет несимметричность решения, в то время как данные [22] предполагают симметричность исходной геометрической модели. По результатам моделирования можно составить итоговую диаграмму распределения поглощенной дозы данной задачи, которая представлена на рис. 11.

Краткие итоги четвертой главы заключаются в следующем.

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

В заключении сформулированы основные результаты работы.

Идея диссертации подсказана научным руководителем и учителем Андросенко Петром Александровичем, посвятившим свою жизнь в науке исследовшшю и применению метода Монте-Карло.

СПИСОК ЦИТИРУЕМОЙ ЛИТЕРАТУРЫ

1. Perkins S., Cullen D., Seltzer S.. Tables and Graphs of Electron-Interaction Cross Sections from 10 eV to 100 GeV Derived from the LLNL Evaluated Electron Data Library (EEDL), Z~ 1-100, UCRL-50400, Vol.31, - Lawrence Livermore National Laboratoiy 1991.

2. Perkins S.T., Cullen D.E., etal. Tables and Graphs of Atomic Subshell and Relaxation Data Derived from the LLNL Evaluated Atomic Data Library (EADL), Z= 1-100, UCRL-50400, Vol. 30, - Lawrence Livermore National Laboratory, 1991

3. ENDF-102, DATA FORMATS AND PROCEDURES FOR THE EVALUATED NUCLEAR DATA FILE, ENDF-6, BNL-NCS-44945-01/04-Rev., Informal Report, Revised April 2001. Written by the Members of the Cross Section Evaluation Working Group Edited by: V. McLane, NATIONAL NUCLEAR DATA CENTER BROOKHAVEN NATIONAL LABORATORY UPTON, N.Y. 11973-5000.

4. MCNP - A General Monte Carlo N-Particle Transport Code, version 4B, March 1997. Edited by Judith F. Briesmeister, Los Alamos National Laboratoiy manual, LA-12625-M, 1997.

5. SalvatF, FemSndes-Vanea J.M., SempauJ.. PENELOPE, a code system for Monte Carlo Simulation of Electron and Photon Transport. // OECD Nuclear Energy Apency Issv-les-Moulineaux, France, 2003.

6. Жуковский M.E., Скачков M.B.. Статистические модели электронной эмиссии Модель «Утолщенныхтраекторий». //Москва, препринт, Институт Прикладной Математики, Российская Академия Наук, 2007.

7. Андросенко П.А., Ефименко Б.А.. Модификация метода Монте-Карло для расчета локальных характеристик потока излучения. И Журнал «Вычислительная математика и математическая физика», 1978, Т.18, №6.-С. 1493-1499.

8. Андросенко П.А.. Применение ПА-метода Монте-Карло для расчета альбедо. // Обнинск, Физико-Энергетический Институт, препринт ФЭИ-975, 1979.

9. Андросенко П.А., Ефименко БА.. ПА-метод Монте-Карло для расчета локальных характеристик потока излучения. // Журнал «Атомная энергия», Т.45, Выпуск 4, 1978. -С. 290-292.

10. Cobut V., Cirioni L., and Patau J.P.. Electron Transport Simulation in the Range 1 keV-4 MeV for the Purpose of High-Resolution Dosimetric Application. // Advanced Monte Carlo for Radiation Physics, Particle Transport Simulation and Applications, Session Electron-Gamma, Proceedings of the Monte Carlo 2000 Conference, Lisbon, October 23-26 2000.

11. Singh J.J. Transmission of 2.43 MeV Electrons Through Thick Silicon Targets. // NASA Technical Note D-5075,1969.

12. Nakat Y. // Jpn. J. Appl. Phys. 2, p. 743,1963.

13. Гусев Н.Г., Машкович В.П., Суворов А.П.. Физические основы защиты от излучений. Издание второе, переработанное и дополненное. // Москва «Атомиздат», 1980г..

14. TorbjOrn Granlund. GNU MP: The GNU Multiple Precision Arithmetic Library. Version 4.1.4,21 September 2004, Swox AB, tege@swox.com.

15. MPFR: The Multiple Precision Floating-Point Reliable Libraiy. Version 2.0.1, The MPFR team, LORIA/INRIA Lorraine, April 2002.

16. Коньков A.B.. Математическое моделирование методом Монте-Карло дозиметрических задач внутреннего облучения в радионуклидной терапии. //Диссертационная работа, Специальность: 05.13.18 - Математическое моделирование, численные методы и комплексы программ, ОГТУАЭ (ИАТЭ), Обнинск, 2008г..

17. Андросенко П.А., Коньков А.В.. Математическое моделирование кровеносной системы для оценки дозовых распределений, создаваемых микросферами, активированным Re-188. И Альманах клинической медицины, Т. XVII, часть 1,. Москва: МОНИКИ. 2008г. - С. 139-143.

18. ХайнДж., БраунеллГ. Радиационная дозиметрия. - Москва: Иностр. лит-ра, 1958,758 с.

19. Stefano Agosteo, Jean-Louis Chartier, Bernd GroBwendt, Gianfranco Gualdrini, Ivan Kodeli, Peter Leuthold, Stdfanie Mdnard, Robert Price, Bernd Siebert, Hamid Tagziria, Rick Tanner, Michel Terrisol, Maria Zankl. QUADOS: "Quality Assurance of Computational Tools for Dosimetry". // Intercomparison on the Usage of Computational Codes in Radiation Dosimetry, Supported by the European Commission, DG XII under Contract FIGD-CT-2000-20062.

20. Robert Alan Price. QUADOS. Problem P2: Endovascular radiotherapy. // Bologna workshop: Intercomparison on the Usage of Computational Codes in Radiation Dosimetry, Italy, July 14-16 2003.

21. Krzysztof Wincel, Barbara Zareba. MCNP Dose Distributions Calculations for 32P Brachytherapy Wire Source QUADOS - P2. Endovascular radiotherapy problem. // Bologna workshop: Intercomparison on the Usage of Computational Codes in Radiation Dosimetry, Italy, July 14-162003.

22. Zilio V.O., Joneja O.P. and ChawlaR.. Problem 2: Dosimetry benchmark of a "P endovascular source. // Bologna woricshop: Intercomparison on the Usage of Computational Codes in Radiation Dosimetry, Italy, July 14-16 2003.

23. Carlos Oliveira and НёНо Yoriyaz. Endovascular Radiotherapy Problem - P2. // Bologna workshop: Intercomparison on the Usage of Computational Codes in Radiation Do-simetiy, Italy, July 14-16 2003.

24. Gualdrini G„ Agosteo S., Minard S., Price R.A., Chartier J.-L., GroQwendtB., Kodeli I., Leuthold G.P., Siebert B.R.L., Tagziria II., Tanner R.J., TerrissolM., Zankl M.. "QUADOS" Intercomparison: A Summary On Photon And Charged Particle Problems. //

21st Century Challenges ill Radiation Protection and Shielding, 1CRS-10, RPS-2004, Madeira, 9-14 May 2004.

СПИСОК ПУБЛИКАЦИЙ ПО ТЕМЕ ДИССЕРТАЦИИ

Статьи в научных рецензируемых журналах

1. Андросенко П.А., Белоусов В.И., МогулянВ.Г.. Модификация метода индивидуальных соударений для моделирования переноса электронов методом Монте-Карло. // Журнал «Вопросы Атомной Науки и Техники», серия: «Физика Ядерных Реакторов», Выпуск 1, Российский Научный Центр «Курчатовский Институт», 2012г., С. 18-23.

2. Андросенко П.А., Белоусов В.И., ЦаринаА.Г.. Исследование методом Монте-Карло воздействия ионизирующего излучения на химический состав тканей человека. II Т.78. Труды регионального конкурса научных проектов в области естественных наук. Выпуск 13. Тематика: Математика, информатика и механика. Калуга: Издательство АНО «Калужский научный центр», 422с., 2008г. - С. 62-73.

3. Андросенко П.А., Белоусов В.И., ЦаринаА.Г.. Моделирование методом Монте-Карло воздействия ионизирующего излучения на химический состав тканей человека. // Журнал «Известия Высших Учебных Заведений. Ядерная Энергетика», Выпуск 3, Раздел: Моделирование Процессов в Объектах Ядерной Энергетики, Обнинский Государственный Технический Университет Атомной Энергетики. 2007г. - С. 84-88.

4. РЛ. Androsenko, V.I. Belousov, A.G. Tsarina. Solution of Transport Problem by Monte-Carlo Technique without Simulating of the Particles Fire Path Length. II Session 6. Monte-Carlo, Probabilistic and Stochastic Transport Methods. Book of Abstracts. The 20th International Conference on Transport Theory (ICTT-20), Obninsk, Russia, July 22-28, 2007.-P. 196-199.

5. Андросенко П.А., Белоусов В.И., Царина А.Г.. Прецизионное решение задач переноса электронов методом Монте-Карло. // Журнал «Известия Высших Учебных Заведений. Ядерная Энергетика», Выпуск 3, Раздел: Моделирование Процессов в Объектах Ядерной Энергетики, Обнинский Государственный Технический Университет Атомной Энергетики. 2007г. - С. 78-83.

6. Андросенко П.А., Белоусов В.И., Коньков АЛ., Царина АГ.. Современный статус комплекса программ BRAND. // Журнал «Вопросы Атомной Науки и Техники», серия: «Физика Ядерных Реакторов», Выпуск 1, Российский Научный Центр «Курчатовский Институт». ISSN 0205-4671, УДК 621.039.5+519.85,2006г. - С. 74-84.

7. Андросенко П.А., Белоусов В.И.. Решение задач переноса электронов методом индивидуальных столкновений. И Математическое моделирование: тез. докл. Региональной студенческой конференции, Обнинск: ИАТЭ, 18-19 мая 2006г. - С. 3-4.

8. Androsenko Р.А., Belousov V.I.. The Monte Carlo sampling of electron transport by individual collision technique. //American Nuclear Society Topical Meeting in Monte Carlo, Chattanooga Marriott and Convention Center, TN, 2005, The Monte Carlo Method: Versatility Unbounded In A Dynamic Computing World Chattanooga, Tennessee, USA, April 1721,2005, on CD-ROM, American Nuclear Society, LaGrange Park, IL, 2005.

9. Androsenko P.A., Belousov V.I.. The Monte Carlo sampling of electron transport by individual collision technique. // Book of Abstracts. International Topical Meeting on Mathematics and Computation, Supercomputing, Reactor Physics and Nuclear and Biological Applications, Palais des Papes, Avignon, France, September 12-15,2005.

10. Андросенко ПЛ., Белоусов В.И.. Прецизионное решение задач переноса электронов методом моделирования индивидуальных столкновений. // Журнал: «Научная сессия МИФИ-2005», Т.5 «Медицинская физика и техника, биофизика. Моделирование фи-

зических процессов в окружающей среде. Охрана окружающей среды и рациональное природопользование. Теоретические проблемы физики», Москва, 2005г. - С. 20-22.

Материалы конференций и тезисы докладов

1. Современный статус комплекса программ BRAND для моделирования переноса нейтронов, фотонов и заряженных частиц. / Всероссийская конференция по вычислительной математике КВМ-2011, Новосибирск, Академгородок! Институт вычислительной математики и математической геофизики СО РАН, 29 июня - 1 июля 2011г

2. Полуаналитический метод Монте-Карло в задачах переноса электронов / Не'й-троника. Ежегодный XVIII семинар «Нейтрошш-физические проблемы атомной энергетики», Обнинск, 30 октября - 2 ноября, 2007г..

3. Solution of Transport Problem by Monte-Carlo Technique without Simulating of the Particles Free Path Length. / ICTT-20. The 20th International Conference on Transport The-

S'july22-28 тГМ' Pr0babilistic md Stochastic Transport Methods, Obninsk, Rus-

4 Electron transport by Monte Carlo individual collision technique. / Workshop on Use of Monte Carlo Techniques for Design and Analysis of Radiation Detector. Symposium topic Modehng and simulation of radiation transport. Coimbra, Portugal, 15-17 September, 2006 '

5. Electron transport by Monte Carlo individual collision technique. / ISRP-10 10th International Symposium on Radiation Physics. Symposium topic: Modeling and simulation of radiation transport. 17-22 September, 2006.

6. Решение задач переноса электронов методом индивидуальных столкновений / III Международная конференция «Математические идеи ПЛ.Чебышева и их приложение к современным проблемам естествознания», Обнинсх, 14-18 мая 2006г

7. Андросенко П.А., Белоусов В.И., Сравнительный анализ результатов расчетов переносов электронов по программам BRAND и PENELOPE. / Нейтроника. Ежегод-

<<НейтР°кно-физические проблемы атомной энергетики», Обнинск

8-10 ноября 2005г.

8 The Monte Carlo sampling of electron transport by individual collision technique / International Topical Meeting on Mathematics and Computation, Supercomputing, Reactor Physics and Nuclear and Biological Applications, Palais des Papes, Avignon, France, Sep-tcmDcr 1 /—i j, 2005.

9 The MonteCario sampling 0f electron transport by individual collision technique. / Monte Carlo 2005. American Nuclear Society Topical Meeting in Monte Carlo, Chattanooga Marriott and Convention Center, Chattanooga, Tennessee, USA, April 17-21 2005

10. Прецизионное решение задач переноса электронов методом моделирования индивидуальных столкновений. / Научная сессия МИФИ, 24-28 января 2005г

11. Моделирование процесса переноса электронов методом индивидуальных столкновений. / Неитроника. Ежегодный XV семинар «Нейтроино-физические проблемы атомной энергетики» 25-29 октября, 2004г.

--- Компьютерная верстка В.И. Белоусов

ЛР№ 020713 от 27.04.1998 ---

Подписано к печати ¿0,03.^, Формат бумаги 60x84/16

Печать ризограф. Бумага МВ Печ л , 5

3аКм№ гх* _Тираж75экз._денадоговорная

Отдел множительной техники ИАТЭ 249035, г. Обнинск, Студгородок, 1

61 12-1/849

Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования ОБНИНСКИЙ ИНСТИТУТ АТОМНОЙ ЭНЕРГЕТИКИ -филиал Национального исследовательского ядерного университета «МИФИ»

На правах рукописи УДК 519.85+614.876+621.039.5

БЕЛОУСОВ ^СЛ

Виктор Иванович г? ^

ПРЕЦИЗИОННЫЕ МЕТОДЫ МОНТЕ-КАРЛО ДЛЯ РАСЧЁТА ТРАНСПОРТА ЭЛЕКТРОНОВ

Специальность: 05.13.18 - «Математическое моделирование, численные

методы и комплексы программ»

ДИССЕРТАЦИЯ на соискание учёной степени кандидата физико-математических наук

Научный руководитель -

доктор физико-математических наук,

профессор Кураченко Ю. А.

Обнинск-2012

Оглавление

Введение...................................................................................................................5

Глава 1. Решение задач переноса электронов...................................................18

1.1. Задачи переноса излучения......................................................................19

1.2. Метод Монте-Карло в задачах переноса излучения.............................20

1.3. Алгоритмическое описание метода Монте-Карло................................22

1.3.1. Особенности моделирования аналогичных программ..................22

1.3.2. Недостатки программ моделирования переноса электронов........23

1.4. Взаимодействие заряженных частиц с веществом................................26

1.4.1. Упругое рассеяние.............................................................................27

1.4.2. Неупругое рассеяние.........................................................................30

1.4.3. Потери энергии заряженных частиц................................................31

1.4.4. Радиационные потери энергии и тормозное излучение................35

1.4.5. Пробег заряженных частиц...............................................................39

1.5. Файлы оценённых ядерных данных........................................................42

1.6. Метод моделирования индивидуальных соударений частиц..............44

1.6.1. Теоретическое описание метода......................................................45

1.6.1.1. Использование интегрального уравнения переноса в методе Монте-Карло...............................................................................46

1.6.1.2. Построение цепи Маркова........................................................47

1.6.1.3. Интегральное уравнение переноса излучения электронов.... 5 0

1.6.2. Моделирование длины свободного пробега...................................52

1.7. Краткие итоги первой главы....................................................................54

Глава 2. Алгоритмическая реализация реакций взаимодействия электронов с

веществом..............................................................................................55

2.1. Использование файлов оценённых ядерных данных............................55

2.2. Основные алгоритмы моделирования взаимодействий электрона.....57

-32.3. Упругое взаимодействие..........................................................................59

2.4. Тормозное излучение...............................................................................63

2.4.1. Моделирование энергии фотона......................................................64

2.4.2. Моделирование энергии электрона.................................................66

2.5. Возбуждение атома...................................................................................68

2.6. Электроионизация атома..........................................................................69

2.7. Релаксация атома......................................................................................71

2.8. Краткие итоги второй главы....................................................................76

Глава 3. Полу аналитический метод моделирования задач переноса

электронов.............................................................................................77

3.1. Применение полуаналитического метода..............................................77

3.2. Реакции взаимодействия электронов с веществом...............................79

3.3. Преобразование уравнения переноса.....................................................80

3.4. Перенос электронов..................................................................................82

3.5. Рекуррентные соотношения.....................................................................86

3.6. Рекуррентные соотношения для барьера заданной толщины..............89

3.7. Рекуррентные соотношения для полупространства..............................91

3.8. Решение проблем вычисления коэффициентов через рекуррентные соотношения..............................................................................................92

3.8.1. Описание библиотеки вМР..............................................................93

3.8.2. Применение библиотеки вМР.........................................................93

3.8.3. Описание расширенной библиотеки МРИ1....................................94

3.8.4. Применение библиотеки МРШ.......................................................95

3.9. Краткие итоги третьей главы...................................................................96

Глава 4. Результаты моделирования задач переноса электронов...................97

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

индивидуальных соударений...................................................................97

4.1.1. Моделирование процесса переноса излучения электронов в кремнии...............................................................................................97

4.1.2. Задача внутрисосудиетой радиотерапии.........................................98

-44.1.2.1. Задача брахитерапии..................................................................98

4.1.2.2. Источник излучения /?-частиц.................................................101

4.1.2.3. Построение модели..................................................................101

4.1.2.4. Оценка дозового поля..............................................................103

4.1.2.5. Результаты моделирования.....................................................103

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

сосудистую систему микросфер....................................................107

4.1.3.1. Аналитическое представление пространственного

188г>

распределения поглощенной энергии для изотопов Ке, 16бНо, 165Dy, 90Y.........................................................................109

4.1.3.2. Эффект самопоглощения излучения в материале микросферы для изотопов 188Re, 16бНо, 165Dy, 90Y.......................................113

4.2. Решение задач переноса электронов полуаналитическим методом.. 117

4.2.1. Моделирование процесса переноса излучения электронов в алюминии.........................................................................................118

4.2.2. Результаты моделирования воздействия облучения на химический состав органической ткани.............................................................119

4.3. Краткие итоги четвёртой главы.............................................................122

Заключение..........................................................................................................123

Список литературы.............................................................................................125

Список рисунков.................................................................................................134

Список таблиц.....................................................................................................138

Приложения.........................................................................................................139

Введение

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

На сегодняшний день математическое моделирование процессов взаимодействия ионизирующего излучения с объектами сложной геометрии и внутренней структуры получило большое распространение. Широкую популярность приобрели вычислительные алгоритмы, основанные на статистическом моделировании методом Монте-Карло [1,2,3,4] процессов переноса и взаимодействия излучения с веществом. Преимущество метода Монте-Карло перед методами, основанными на численном решении кинетического уравнения, определяется его удобством и приспособленностью к решению сложных граничных задач в многокомпонентных средах [5, 6, 7].

Исследование возможностей использования метода Монте-Карло для решения задач переноса электронов определило направленность данной работы.

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

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

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

Повышение точности выполняемых расчётов электронно-физических характеристик в большой степени помогает решению рассматриваемых задач. Требования повышения точности выполняемых расчётов диктуют, в свою очередь, необходимость использования самой современной информации о взаимодействии излучений с веществом, которая содержится, как правило, в файлах оценённых ядерных данных (например, зарубежные библиотеки ЕЕОЬ [8, 9], ЕАБЬ [10, 11], ЕРБЬ [12, 13], ШЖ>Ъ-3 [14], РЕЖ>Ь-2 [15, 16, 17], отечественная В1ЮЖ>[18]). Обеспечение необходимой точности при решении уравнения переноса излучения возможно, как правило, лишь при подробном описании реальной трёхмерной геометрии исследуемого объекта и при детальном учёте информации о взаимодействии излучения с веществом. Все вышесказанное наиболее корректно может быть выполнено при использовании программных

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

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

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

Первоначальные методы моделирования процесса переноса электронов [19] в основном состояли в использовании табличной информации о прямолинейном пробеге электронов [20]. Дальнейшее применение метода Монте-Карло наряду с дополнением табличной информации сечениями электронно-атомных взаимодействий привело к созданию новых алгоритмов [1,5] переноса электронов в веществе. Полученные алгоритмы были использованы в разработках новых программ моделирования процесса переноса электронов. Например, в программном комплексе МС№ [1], где используется модель укрупнённых историй.

Модификацией данной модели является модель утолщённых траекторий [21], эффективно применяемая для задач глубокого проникновения электронов. Существует комбинированная модель укрупнённых соударений с детальным описанием процесса взаимодействий электрона с атомами среды, которая реализована в виде программного комплекса PENELOPE [2].

В рамках программных комплексов, реализующих использование метода Монте-Карло для решения задач переноса излучений, работу по моделированию процессов взаимодействия частиц с веществом выполняют подпрограммы так называемого константного модуля, которые извлекают необходимые данные из соответствующих библиотек оценённых данных. Данный модуль является одной из самых трудоёмких частей монте-карловского комплекса как с точки зрения временных затрат при расчёте, так и в смысле математического моделирования, физического обоснования, алгоритмизации и программной реализации. Специфика методов с константным модулем в монте-карловских программах позволяет использовать имеющуюся информацию о взаимодействии излучения с веществом практически без упрощений, вплоть до прямого извлечения из файлов оценённых данных. Использование современных библиотек оценённых данных по электронно-атомным взаимодействиям EEDL дало возможность применять методы, разработанные на основе метода Монте-Карло, с дополнительным увеличением точности этих методов.

Цель работы состояла в дальнейшем развитии комплекса программ BRAND [3], как прецизионного инструмента для численной имитации оценочных экспериментов. В рамках данного программного комплекса, используя полученный модуль подготовки константной информации по файлам библиотек оценённых данных EEDL и EADL, была произведена разработка и программная реализация алгоритмов моделирования процесса переноса электронов.

Метод Монте-Карло является единственным методом, позволяющим точно моделировать процесс переноса излучения. Перенос электронов

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

Использование расчётной программы моделирования переноса излучения электронов связано с необходимостью решить проблему выбора констант, определяющих точность решения задачи переноса. В нашем случае эта проблема решается с помощью библиотек файлов оценённых данных EEDL [8] и EADL [10]. При использовании других численных методов моделирования переноса излучения существует вероятность потерять некоторые важные данные о характере траектории электронов в среде и получить результат без их учёта. Предложенный метод моделирования индивидуальных соударений [24,25,26,27] позволяет обойтись без предварительных предположений, а также дополнительных осреднений потока частиц [1,2,21], и начать поиск решения задачи, применяя классическую модель метода Монте-Карло, без предварительных исследований, используя «напрямую» информацию из файлов оценённых данных формата ENDF-6 [23] библиотек EEDL [8] и EADL [10], что позволяет ограничить точность получаемого численного решения только данными константного модуля.

В рамках комплекса моделирующих программ BRAND была произведена разработка и программная реализация алгоритмов моделирования процессов рассеяния электронов в процессе монте-карловского расчёта по информации файла 26 формата ENDF-6 оценённых данных библиотеки EEDL «напрямую», без внесения каких бы то ни было приближений и упрощений по реакциям взаимодействия электрона с атомом вещества. Была также произведена дополнительная разработка и программная реализация алгоритмов моделирования процесса релаксации атома «напрямую» по информации из файла 28 формата ENDF-6 библиотеки EADL после воспроизведения процесса электроионизации атома.

-10В результате был создан уникальный константный модуль программного комплекса BRAND, в котором удалось объединить всю необходимую информацию о взаимодействиях электронов с атомами вещества. Проверка полученного константного модуля и алгоритмов моделирования взаимодействий электрона с материалом были сделаны благодаря применению полуаналитического метода (ПА-метод) [28, 29, 30], который позволяет имитировать энергопотери частицы в веществе без моделирования длины свободного пробега, что существенно уменьшает погрешность вычислений и позволяет получать фактически аналитическое решение задачи переноса.

Научная новизна настоящей работы заключается в следующем:

1. Разработана трёхмерная математическая модель имитации задач переноса электронов - метод индивидуальных соударений, который является разновидностью классического моделирования методом Монте-Карло.

2. Созданы и программно реализованы алгоритмы моделирования процессов рассеяния электронов по информации файлов оценённых ядерных данных библиотек EEDL и EADL формата ENDF-6.

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

4. Разработана модификация версии программного комплекса BRAND для решения задач переноса электронов с помощью полуаналитического