Обоснование DXTRAN-модификации метода Монте-Карло на основе соотношений взаимности для различающихся систем Текст научной статьи по специальности «Математика»

УДК 519.245

Е.А. Цветков1,2; В.В. Шаховский2

1 Московский физико-технический институт (государственный университет)

2 Центральный научно-исследовательский институт химии и механики

Обоснование DXTRAN-модификации метода Монте-Карло на основе соотношений взаимности для различающихся систем

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

Ключевые слова: перенос ионизирующего излучения, метод Монте-Карло, весовые методы, соотношение взаимности, DXTRAN, MCNP.

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

Одним из методов расчёта переноса ионизирующего излучения является метод Монте-Карло. В случае малого по своим геометрическим размерам детектора он требует розыгрыша очень большого числа траекторий. Причиной этого является низкая вероятность попадания частиц в область детектора. В наиболее простом случае точечного детектора можно применить метод локальной оценки, в котором при розыгрыше столкновений частица переносится сразу в точку, где находится детектор, причём её вес рассчитывается с учётом вероятности рассеяться в направлении этой точки и с учётом оптического пути при прохождении через слой вещества к этой точке. Для распределенной системы детекторов в работах [1] и [2] предложен метод решения подобных задач через перевязку прямого и сопряженного решений уравнения переноса на поверхности, окружающей интересующую нас область. На этой поверхности случайным образом выбирается большое количество точек. Из

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

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

Энергия DXTRAN-частицы определяется после розыгрыша её направления.

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

В документации [8] к программе MCNP метод DXTRAN обосновывается подведением баланса весов. В то же время нам представляется возможным обосновать указанную модификацию на основе общих принципов, привлекая для этого соотношения взаимности для различающихся систем, представленные в работе [3] для нулевых краевых условий и обобщённые на краевые условия общего вида в работах [5, 6]. Для расчёта методом Монте-Карло функционалов от потока внутри локализованной области (в нашем случае внутри DXTRAN-сферы) в работе [9] был предложен алгоритм, в котором прямыми блужданиями моделировалось решение Ф* (ж) сопряженной задачи.

Для лучшего понимания рассматриваемого метода повторим сначала обоснование метода подведением баланса весов, приведённое в [8]. Это обоснование не может считаться строгим, но помогает понять метод интуитивно.

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

Пусть w0 — вес частицы до столкновения. Пусть р1 — вероятность того, что участок траектории частицы после разыгрываемого столкновения и до следующего столкновения не пересечет DXTRAN-сферу, а р2 — вероятность того, что этот участок пересечет сферу. При соударении частица расщепляется на две частицы — одна с весом w0p1, вторая с весом w0p2. На этапе расщепления сумма весов двух частиц остаётся прежней: w0p1 + w0p2 = w0.

Величины р2 и р1 определяются по формулам

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

Новая частица с весом w0p1 рассеивается обычным способом. В случае если она пересекает сферу до следующего столкновения, то эта частица уничтожается. Можно сказать, что на этом этапе происходит игра в русскую рулетку с вероятностью выживания р1. В случае выживания в русской рулетке вес частицы должен быть увеличен в 1/р1 раз и составит w0p1/p1 = w0, то есть выжившая частица продолжает двигаться с весом, равным весу исходной частицы до столкновения. Усреднённый вес этой частицы равен

и)1 = w0p1.

Заметим, что величина p1 неизвестна.

Новая частица с весом w0p2, называемая DXTRAN-частицей, рассеивается в направлении на сферу и переносится в соответствующую этому направлению точку на поверхности сферы. Рассеяние разыгрывается с вспомогательной плотностью вероятности Рагь(П), причём на всех направлениях, не указывающих на сферу, Рагь (П) = 0, что эквивалентно

Рагь(П)в,П = 1.

Далее применяется обычная схема для расчёта веса частицы, рассеянной со вспомогательной плотностью вероятности. А именно, вес частицы умножается на величину

Р(П)

1^2

1)1 = 1 - 1)2,

(1)

Р(П\П е Б)

Рагь(П) РогЬ(П) /Р(П)сЮ’

равную отношению условной плотности вероятности Р(П|П € Б) рассеяться в направлении П при условии, что рассеяние произошло в направлении на сферу к вспомогательной плотности вероятности. При

переносе на сферу вес частицы умножается ещё на в-т(п). Окончательно получаем, что вес DXTRAN-частицы равен

W2 = Wо

Р (П)в-т(п) йПх

X

Р (П)в-т(п)

РагЬ (П) / Р(П)в-т(П) йП

УО о-

Р(П)е~т(-п) Рагь(П)

(2)

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

Определим среднее значение веса DXTRAN-частицы:

и}2

Ргке-т{П)РагЬ(П)(1П.

РагЬ (П)

С учётом (1) получаем

0)2 = WоP2.

Из вышеприведённых выкладок следует, что в методе DXTRAN сумма весов частиц в среднем сохраняется:

У01 + й)2 = Wо,

причём не возникает необходимости вычислять интеграл в (1).

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

Итак, рассмотрим две задачи а и Ь, определённые в области и с границей Г.

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

Рис. 1. Геометрия задач а и Ь

В работе [6] показано, что справедливо выражение

(К,фь)г = Ра$д — Рь8р — (Ра — Рь)8М +

+ (Фа,Фа)Г 8Я + (Фа,Фь — Фа^ 8М. (3)

В нем символами 8д, 8Р, 8м обозначены индикаторы попадания подобластей ид, ир и им соответственно в область и', величины Ра и Рь являются откликами детекторов в задачах а и Ь соответственно и определяются формулами

Ра

Рь

Фа(г,ш,Е ^(г^Е )йЕйшйг,

Фь(г,ш,Е ^(г^Е )йЕйийг,

и (4п) Е

в которых интегрирование ведётся по всей области и, по всем направлениям в пределах телесного угла 4п и по всем энергиям в пределах от 0 до максимально возможной энергии в данной задаче. Выражение а(х,и,Е),д(х,и,Е))г обозначает интеграл по энергии в тех же пределах, по направлениям в тех же пределах и по поверхности Г:

а (х,^Е )д (х,^Е ))г = а (х,и,Е )д(х,и,Е )(п,ш)йЕйшйх,

Г (4п) Е

где п — внешняя нормаль к поверхности.

Пусть границей Г' является DXTRAN-сфера. Пусть задача Ь получается из задачи а введением абсолютного поглотителя в области и'. Индикаторы в этом случае принимают значения: 8„

0, 8Р

1.

8м = 1- Тогда выражение (3) упрощается: <Ф;,ФЬ)Г, = -Fa + {Ф1,Ф„ - Фа)г . (4)

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

Фь- = РЬ + Rb Фь+;

Ф* = D* + R* Ф

Фа+ D a + 1ьаФа-,

(5)

Фа- — Ра + КаФЬ+ ■

Оператор отражения и сопряженный оператор отражения связаны формулами [10]:

RФb ,Фа>г_ +

Фь,я *ф:

г+

КФ,ФА

+ (Ф,Д * Ф* г^ \ /г+

0,

0. (6)

Перепишем второй член правой части выражения (4) с учётом (5) и (6) в виде

(Фа.Фь - Фа>г —

— Р.фь - фа>г+ + (фа.Рь - Ра>Т- ■

Правая часть этого выражения равна нулю, если потребовать выполнения равенств

' ра — о.

Рь — Ра.

Первое условие означает, что за пределами области и детекторы отсутствуют, второе условие требует, чтобы в задачах а и Ь источники за пределами области и создавали одинаковые потоки на границе Г. Эти равенства в случае решения задачи методом DXTRAN выполняются. В этом случае формула (4) приобретает вид

<Ф'аФь}Р = -Fa-

(7)

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

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

источником, до столкновения в точках хг. Буквами уг^ — фазовые координаты DXTRAN-частиц до столкновения в точках уг,^. Буквами и ^ — направления движения этих частиц до столкновения, Ег и Ег^ — энергии частиц до столкновения. Можно записать, что хг = (хг,шг,Ег) и уг,з = (уг,з М,з Ег, з). Чтобы сократить объём выкладок везде, где это не приведёт к путанице, будем опускать дельта-функции вида 8(Е — Е') и интегрирование по энергии.

Плотность излучения Фь(ж) в задаче Ь удовлетворяет уравнению

Фь(х) — КьФь(х) + Фоь (ж), (8)

где Кь — интегральный оператор переноса и Ф0ь(ж) — плотность нерассеянного излучения. Функция ценности Фа (ж) в задаче а удовлетворяет уравнению

Фа(ж) — ка Фа(ж)+Ф0а м. (9)

где Ф0а(ж) — нулевой член разложения по столкновениям функции ценности, Ка — сопряженный оператор переноса. Решение каждого из этих уравнений можно представить в виде ряда Неймана, который мы должны получить ниже.

Усредним по всем возможным траекториям случайную величину, равную отклику детектора на DXTRAN-частицы, порожденные одной и той же частицей в задаче Ь. Эта величина равна сумме откликов детектора ф(уі,з) на акт взаимодействия в точках уі^, і > 1, умноженных на соответствующий вес DXTRAN-частицы:

п Г>( \ к(п)

j=i

где т(хг ^ уг,0) — оптический путь из точки хг в точку уг,0, с энергией частицы, соответствующей энергетической компоненте фазовой переменной уг,0, Р(хг ^ уг,0) -дифференциальная по телесному углу плотность вероятности того, что частица рассеется в точке хг в направлении на точку уг,0, соответствующая реальному физическому процессу, а Рагь(хг ^ уг,0) — плотность вероятности, с которой моделируется рассеяние. Верхние пределы суммирования являются случайными величинами. Чтобы выполнить суммирование корректно, следуя [11], дополним сумму до бесконечной, введя случайную величину:

Д.

і,з

0,

иначе.

Эта величина равна 1, если в задаче b частица достигла точки Xi, а в задаче a DXT RAN -частица достигла точки yij . В противном случае эта величина равна 0. Теперь выражение для £ можно записать в виде двух бесконечных сумм:

t = ул ул P(Xi ->■ У г,о) с-т{хі-+шл) ^ л Рагь(%і —> У г, о)

г=0 1 = \ 4 ^ /

X

хф<УМ )Дч ■

(10)

Величина

Р(хг Уі,о) с-т(хі->Уі.п)

Рагь(%і ~У і,о)

неотрицательна. Величина Д^- тоже неотрицательна по определению. Величина же Ф(Уі,і) является средним откликом детектора на рассеяние частицы в точке уі ^. Пока примем эту величину тоже неотрицательной. Тогда если ряд (10), все члены которого неотрицательны, сходится, то он сходится абсолютно и справедливо равенство

ОО ОО

мц уу М(^Щ х

j==l Parb (xi ^ Уі.0)

Xe-m(xi—yi Myhj )Ді; )■

Далее выделим из суммы член с i = 0:

М£ = м( ~* ^°’°) с~т(хо^Уо,о) у

\Рагь(х0 У0,0)

X Ф(у0.j )Д0 . А +

j=1 X <х

+уум1 х

^ ^ ' Parb(Xi ^ Уі.о)

i=1 j=1

e

-m(xi —>yi, о

]Ф(УіЛ )Д

j/^i. j І '

(11)

Вычислим математическое ожидание произвольного члена ряда в (11) при г > 1. Для этого используем формулу полного математического ожидания:

M

)ф(Уі^ )Ді,з

Р(Хі У і, о) m(.Tj—>'^¿.0

Parb(x'i ^ Уі,о)

Р(*г Уі,о)

Х^т’тА-"ШЛРагЬ(Хг ^ Уі,0)

Xe-m(xi—yi, о)ф(уі j )X

Mx

X

xM(Ді. j \хоХі...ХіУі.оУі. 1 ...уі. j)

(12)

Обозначим д(х') — вероятность поглощения в точке х', г(х' ^ х) — условная плотность вероятности того, что частица,

ется в направлении точки х и достигнет её при условии, что рассеяние произошло. Обозначим через к(х' ^ х) произведение г(х' ^ х)(1 — д(х')), имеющее смысл вероятности того, что частица, испытав столкновение в точке х', выживет, рассеется в направлении точки х и достигнет ее. Эти величины будем снабжать нижним индексом а или Ь в зависимости от того, к какой задаче они относятся. Величина М(Аг,з 1х0х1...хгуг,0уг, 1...у%,з) имеет смысл вероятности того, что частица, испущенная в задаче Ь в точке х0, не будет поглощена ни в одной из последующих точек х1, ..., хг-1, а частица, полученная переносом DXTRAN-частицы из точки хг в точку уг 0, не будет поглощена ни в одной из точек уг, 1, ..., уг, з-1 в задаче а и достигнет точки уг з. Очевидно, что эта величина равна произведению вероятностей выживания частицы при взаимодействии в каждой из этих точек:

М(Дг,з 1х0х1...хгуг,0уг, 1...уг,з) =

і1

j-1

Д (1 - дь(хі)^ П (1 - da^im)).

і=і

m=1

испытав столкновение в точке х . рассе

Для вычисления М£ найдём теперь плотность вероятности реализации цепочки х0х1...хгуг,0уг,1...уг,з. Пусть рождение частицы в задаче Ь разыгрывается с плотностью вероятности Б0(х). Тогда плотность вероятности реализации части цепочки, относящейся к задаче Ь, равна

Р(х0 * Х\

■■■

£—т(хо^хі) Хі) = БоІХо)-----------“X

|Хі - Х0\

хё(ш0 -

Хі - Х0 |*і - *о|

іі

) Д Гь(хі * Хі+і)■

і=і

Условная плотность вероятности перехода из точки хг в точку уг,0 при условии, что в точке хг произошло рассеяние, равна

Рагь(хг ^ уг,0).

Вероятность того, что в точке хг произошло рассеяние, равна

Рз(Хі)

У*(хі)

Рг(Хі)

Плотность вероятности DXTRAN-частице попасть из точки хг на сферу в окрестность точки уг)0 при разыгрывании Рагь(хг ^ уг,0) выражается формулой

дР(хг ^ уг,0) _ , , дП

РагЪ\рС% * Уi,o)~

%

4,0

г,0

где

дп

дуі,0

телесный угол, под которым

видна элементарная площадка на поверхности сферы, содержащая в себе точку уг,0, из точки хг. Величину Р(хг ^ уг$) представим в виде

Р (Хі * у і а)

(?з(хі ~» Уі,о)

Єв(Хі)

где а3(хг ^ уг,0) — дифференциальное по телесному углу сечение рассеяния частицы в точке хг в направлении на точку уг,0 с энергией после рассеяния Ег^0, а3(хг) — полное сечение рассеяния в точке хг.

Плотность вероятности реализации части цепочки, относящейся к задаче а, рав-

на

Р(уі,0 * уі,1

уі,3 ) —

е—т(уі,о ^Уі,і)

\УгЛ ~ Уі,012

X

хё І Ші,0 -

У г ,1 ~ У г,О \'Уг,1 ~ Уі,0І

з—і

т=1

Здесь множитель

1

Іуі,1 - уі,0\

ё I 0 —

У і,і - У І,о \Уі,і - Уі,о\

учитывает тот факт, что координата уі,і лежит на продолжении луча Хіуі,0. Теперь можно записать в виде

Щі,з —

dx0dx1 ...¿хіх

хохі ...Хі

X

dyi,оdyi,l...dyi,j х

уі, 0уі,1..-уі,з

х8о(хо)-

0—т(хо^хі)

ІХі - Х0І

-5{іОо —

Хі - Х0

\хі - ж0|

)х

іі

X

II о : о - х1+1у^^^е-т^-о)х

&8 (Хі )

і=і

X

а8(хі) дП е~т{:Уі^Уіл) О і (хі) дуі, 0 ІУ і, 1 — У і,о І'

у і,і - уі,0

X

хё І ші,0 -

\уі,і - уі,0\

X

з—і

х Д Га (уі,т * уі,т+і)Ф(уі,3)х

т=і

і— і

з—і

Х П (1 — дь(х^ П (1 — да(уг,т)) .

1=1 т=1

Здесь интегрирование по уг,0 производится таким образом, что пространственная часть фазовой переменной уг,0 является точкой DXTRAN-сферы. Объединяя произведения с одинаковым индексом и учитывая, что

р-т(х'^х,Е) х _ х

80(х',и,Е)—г~,------¡—¿(^ ~ т~,-------Мх' =

| х' — х| | х' — х|

У

= Ф0ь(х,ш,Е), проинтегрируем по переменным х0, х1, ..., х—1 и уг,з, уг,з-1, ..., у*д. С учётом

г- 2

dx1 ..Ахі— і Д кь(Хі * х1+і)х

хі ...хі— 1

і=і

хФ0ь(хі) — КЬ-іФь(Хі)

dyi,l...dy

?—т(уі,о^Уі,і)

і,з

уі,і...уі,3

хё ( Ші,0 -

Іуі,і - уі,0\ у і,і - Уі,0

■X

Теперь вернемся к первому члену формулы (11). Аналогично выводу (12) полу-

чаем

\Уі,і - Уі,0\

х

з—і

X П ка (Уіт * Уі,т+іШУі,3) — Ка іФ+

р(хо —» Уо.о) г-т(Х0^У0 0)

Р(Х0 * У0,0)

м , ,)Д0і .

УатЬ (х0 * У0,0)

— М

х0уі,0уі,і...уі,з

РатЬ(х0 * У0,0)

X

т=і

запишем результат в виде

ыс

і,3

dxi

dyi,оdyi,lКЬ іФ0ь(хі)х

хе-т(хі^уі:о)СГз(Хі У г,о) дії

хе т(хо~"У0,0)ф(Уі,з)Ы(А0,з ІХ0У0,0У0,і...У0,3 )

Математическое ожидание величины Д0,3-запишется в виде

з-і

<Ті(Хі)’’ Оуіак° ф0»(«м»)- к0,!/0.0У0,1---У0,у) II I - 9а(Уі.т))

(13)

Величину запишем в виде

дії _ дії — (п,Шуі0)dГ/

дуг,0 д Г' 1хг — уг0

где \хг — уг,0\ — расстояние между пространственными компонентами точек хг и уг,0, п — внешняя нормаль к DXTRAN-сфере Г', а 0 — компонента, соответствующая направлению скорости для фазовой точки уг 0. Знак минус учитывает то, что нормаль п направлена наружу сферы, а 0 направлено внутрь сферы. Величина же должна быть положительной. Формула (13) примет вид

ыс

і,3

dxi

dyi,оdyi,lКгь іФ0ь(хі)х

х е-т{х^у1'0) (Та(хг ~У г, о) ^ сп(*ч)

— (п,шУ. ПЫГ' „

\хг уг,0\

Интегрируя по йхг, получаем ещё одно ядро Кь. Окончательно можно записать:

ыс

і,3

КгьФ0ь(г,ш,Е )К*з—іх

(4п) Е Г'

хФ*а(г,ш,Е) (п,ш) йшйЕйГ =

= — ( КФи,К3-1Ф*^г.

Напомним, что в этой формуле индекс г больше нуля.

т=і

Величину Р(х0 ^ у0,0) представим в виде

е~т(хо ——Уо,о)

Р(хо Уо,о) = *5(^0 —> Уо,о)п------------------гх

\у0,0 — х0\

Х( У°-° ~ х° ^ хё(и)о - -------------

\у0,0 — х0\

где Б(х0 ^ у0,0) — плотность вероятности того, что источник испустит в точке х0 частицу в направлении у0 0. Проведя выкладки, аналогичные случаю г > 0, получаем

ыс

0,3

Ф0ь(г,ш,Е )К*3—іх

(4п) Е Г'

хФ0а(г,ш,Е) (п,ш) dшdEdГ

= — ^,КЗ 1Ф0а) .

Теперь математическое ожидание вклада в детектор можно записать в виде ряда:

М£ = — и и (К'ьФ0ь,К'а’-'Ф0,а)

г=0 з=1 Г

Учитывая, что решения уравнений (8) и (9) можно представить в виде ряда Неймана, окончательно запишем:

М£ = — (Ф;,Фь)р ,

что с учётом (7) даёт

М£ = Ра.

Теперь перейдем к случаю, когда величина ф(уг,з) знакопеременная. Представим её в виде разности:

Ф(уг,з) = Ф+(уг,з) — ф-(уг,з),

2

где

ф+ (уг,з) = тах(0,ф(уг,з)),

ф-(уг,з) = тах(0, — ф(уг,з)).

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

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

1. Разыграть точку источника, которая испустит очередную частицу. Начальный вес новой частицы т0 равен 1.

2. В соответствии со вспомогательной плотностью вероятности Рагь(П) выбрать направление испускания DXTRAN-частицы и разыграть её энергию. Перенести её вдоль этого направления на DXTRAN-сферу. Определить по формуле (2) её вес.

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

4. Если участок траектории до точки следующего столкновения пересекает DXTRAN-сферу, то перейти к п. 8. Иначе перенести частицу в точку следующего столкновения.

5. В точке, куда была перенесена обычная частица, разыграть направление рассеяния DXTRAN-частицы, затем разыграть её энергию в соответствии с углом рассеяния. Определить по формуле (2) вес DXTRAN-частицы. Перенести её на DXTRAN-сферу.

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

7. Перейти к п. 4.

8. Произвести розыгрыш траекторий всех созданных DXTRAN-частиц обычным способом, как если бы они были испущены поверхностным источником на DXTRAN-сфере.

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

Аналогично формуле (7) можно вывести ещё одну формулу, связывающую перевязку решений двух задач и отклик детектора одной из задач. Пусть задача а получается из исходной задачи введением абсолютного поглотителя вне DXRTAN-сферы, а задача Ь совпадает с исходной. Из (3), (5) и (6) в случае 6д = 0, 8Р = 1 и $ы = 0 получаем равенство

<ф; ,Фь>г, = —Рь,

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

Литература

1. Золотухин В.Г., Ксенофонтов А.И., Гнутиков А.П. Решение задач переноса гамма-излучения методом Монте-Карло с помощью теории возмущений // В сб. статей Вопросы дозиметрии и защиты от излучений. — М.: Атомиздат, 1981. Вып. 20. — С. 7-10.

2. Гнутиков А.П., Ермолаев И.В., Золотухин В.Г. [и др.].]. Использование возможностей монте-карловских программ расчёта возмущений от локальных неодно-

родностей // Тезисы доклада на VII Всесоюзном совещании «Методы Монте-Карло в вычислительной математике и математической физике». — Новосибирск, 1985. — С. 252-254.

3. Петров Э.Е. О применении билинейных функционалов для расчёта эффектов возмущения среды в задачах переноса излучения / Препринт ФЭИ-2026. — Обнинск: ФЭИ, 1989.

4. Коробейников В.В., Усанов В.И. Методы сопряжения в задачах переноса излучений // Физика и техника ядерных реакторов. В. 39. — М.: Энергоатомиздат, 1994.

5. Шаховский В. В. Соотношение взаимности для различающихся систем и эквивалентные краевые условия в методах сопряжения задач переноса излучений // Тезисы доклада на VII Российской научной конференции по защите от ионизирующих излучений ядерно-технических установок. — Обнинск: ГНЦ РФ ФЭИ, 1998. — С. 26-28.

6. Шаховский В.В. Формулы метода сопряжений для задач переноса излучений с краевыми условиями общего вида // Тезисы доклада на VIII Российской научной конференции «Радиационная защита и радиационная безопасность в ядерных

технологиях». — Обнинск: ГНЦ РФ ФЭИ, 2002. — С. 21-23.

7. Los Alamos Radiation Transport Group X-16. MCNP-A. Generalized Monte-Carlo Code for Neutron and Photon Transport. — New Mexico: Los Alamos Scientific Laboratory, 1981.

8. RSICC computer code collection. MCNP4C. Monte Carlo N-Particle Transport Code System. — New Mexico: Los Alamos National Laboratory, 2000.

9. Золотухин В.Г., Ксенофонтов А.И., Гнутиков А.П. Новый подход к оценке линейных функционалов и решению задач теории возмущений методом Монте-Карло. Методы Монте-Карло в вычислительной математике и математической физике // Труды VI Всесоюзного совещания. Ч. 1. — Новосибирск, 1979. -С. 163-172.

10. Гермогенова Т.А. Локальные свойства решений уравнения переноса. — М.: Наука, 1986.

11. Ермаков С.М., Михайлов Г.А. Статистическое моделирование. — М.: Наука, 1982.

Поступила в редакцию 04.03.2008.