Алгоритмы метода Монте–Карло для исследования временной асимптотики потока частиц с размножением в случайной среде Текст научной статьи по специальности «Математика»

АПВПМ-2019

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

Г, А. Михайлов1,2, Г, 3, Лотова1,2

1 Институт, вычислительной математики и математической геофизики СО РАН, 630090, Новосибирск 2 Новосибирский государственный университет, 630090, Новосибирск

УДК 519.676

Б01: 10.24411/9999-016А-2019-10055

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

Ключевые слова: метод Монте-Карло, весовые модификации, теория переноса, временная постоянная размножения частиц.

Вводная информация

Рассматривается односкоростной процесс переноса частиц с изотропным рассеянием. В качестве математической модели процесса используется цепь Маркова [хп}, п = 0,1,..., М, соответствующих "столкновений" в фазовом пространстве X = К х V х Т координат, скоростей (с фиксированным значением V модуля) и времени, т.е. хп = (г„, чп,гп), причем \п = у(гп - г„_1)/|г„ - г„_1|, ¿„ = ¿„_1 + |г„_1 - гп|/г>. Эта модель определяется плотностью $(х) точки х0 и плотностыо к(х', х) перехода го состояния ж' в х, причем

/ц^ь = < 1 -4 {><и)

X

и тем самым среднее число переходов конечно, например, в ограниченной системе [1]. Заданы ае(х) и af (х) — коэффициенты рассеяния и деления, среднее число V частиц в точке деления, ас(х) — коэффициент поглощения. Полный коэффициент ослабления потока частиц а(х) = ае(х) + о^(х) + ас(х) (см. [1]).

Методы Монте-Карло используют для оценки линейных функционалов вида ^ = к € Ь^, где

ф(х) — плотность столкновений, удовлетворяющая уравнению = К+ /, в котором К — интегральный оператор с ядром к(х', х). При выполнении условия (1.1) имеем ||_К"< 1 и, следовательно, спектральный радиус оператора р(К) < 1 (см.[1]). Для построения весовой модификации алгоритма используют цепь Маркова с начальной плотностью /0(х) и плотностью перехода р(х',х), вспомогательные веса

N

Яо = ](хо)//о(хо), Qn = Qn-lk(xn-l,xn)/p(xn-l,xn) и "оценку то столкновениям" £ = ^ Япк(хп). Если

п=0

выполняются "условия несмещенности" [1], то Е£ = (^, к), Если щи этом р(Кр) < 1, где Кр — оператор с ядром к2(х',х)^р(х',х),и /2//о € Ь1(Х), то < [1].

Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (коды проектов 18-01-00599, 17-01-00823, 18-01-00356).

!ЯВ.\ 978-5-901548-42-4

Пусть о(х; го, Уо) — плотность столкновений (по аргументу х) от одного столкновения в точке (го, Уо, 0), т. е. для /(х) = 6(г - го)б(у - Уо)6(Ь). Функционал

Д У

можно представить [2] в виде

3(1) = ] ] J ?(го, Уо,* - Т)Р(го, Уо,т)(1 г^Уо(1т,

Д V о

где

Р(го, Уо^) = ! J Ро(г, У, Ц го, Уо)^(г, у)(1г (1у.

Д V

Предполагается, что /(г, у, -) = 0, и вследствие этого Р(г, у, -) = 0, при Ь > 0. Символом /(т) далее обозначается то-кратная производная от функции ] по Лемма 1. [2] Пусть точка (го, уо) распределена для ¿о = 0 с плотностью /о(г, у),

Ц (т)(г, у, Ь)//о (г, у) | <С< +то, т = 0, !,...,п, р(Кр) < 1 и Р(х) < С < +то. Тогда выполняется соотношение 7(т)(г) = Ее(т), где

N

е(т) = $пЧгп, Уп)/(т)(го, Уо^ - и)//о(го, Уо), до = 1,

п=о

причем < +то, т = 0,1,... п.

1 Оценка параметра экспоненциальной асимптотики по времени

Известно, что при выполнении довольно общих условий в случае источника, локализованного в точке (го, Уо, 0), имеет место асимптотическое при Ь ^ то соотношение [3,4]:

^ (г, у, г) - С (г, у)еЛ, С (г, у) <Со < +то, (2.1)

где Л — ведущее характернстпческое число соответствующего однородного стационарного уравнения переноса с заменой ас ^ ас + Л^и. Эти условия, в частности, имеют место для одпоскоростпого процесса переноса частиц в ограниченной среде с достаточно быстро убывающей по времени плотностью источника. Теорема 1. Если выполняются соотношения

у7(п)(г, У,ф-ЛЧг< +то, п = 0,1, (2.2)

соотношение (2.1) и условия леммы 1, то

3 '(*)

^ ^и Ь ^ +то. (2.3)

J (г)

Доказательство. Прямое интегрирование с учетом соотношений (2.1), (2.2) дает равенство

з (г) = Сем [1 + е(г)], = с\ем [1 + £1(г)],

причем е(1), £]_(£) ^ 0 для Ь ^ +то. Интегрируя функцию 7'(¿) в пределах (т, т ^ то для Л < 0

получаем Л'(Ь) = СЛеЛ[1 + £1(£)], т. е С1 = ЛС. В случае Л > 0 тот же результат получается путем введения дополнительного поглощения с коэффициентом тЛг. Это завершает доказательство теоремы 1.

2 Оценка величин ЕЛ и Б Л для случайной среды

Определяемую леммой 1 и теоремой 1 оценку метода Монте-Карло можно рандомизировать с целью изучения флуктуаций значения А для процесса переноса частиц в случайной среде [2]. Будем полагать а(г) — случайное поле, причем отношения а^а, а^а, а также пндпкатрнсы рассеяния и деления фиксированы. При этом J(Ъ) = J(Ъ,а), J'(t) = ,а) и А(а) к (Ь,а) соответственно теореме 1. В дальнейшем

аргумент I будет опускаться, т.е. J(Ь,а) ^ J(а).

Практически важными являются величины ЕА(а) и БА(а). Логически простейший (прямой) сппсоб их оценки методом Монте-Карло состоит в реализации достаточно точных оценок величины J'(t,а)/,](1 ,а) а

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

Е ( ) к А„ -Е J (а)

' ( — 1)8 J' ИЕ^Т ^ (а)

8=0

где J0 к EJ(а).

Простейшая (элементарная) несмещенная рандомизированная оценка величины Ап на основе леммы 1 строится путем реализации п +1 условно независимых траекторий частиц для выбранного значения а : Ап — ЕАп, где

п (— 1)п я

Ап — И'(Пп+и а) £ П (№к,а) - .7о). (3.1)

8=0 и0 к=1

ментов из группы 01:... О 0^ + ^ ...0п+1 после подстановки 0^ ^ 0п+1 в (3.1), в связи с тем, что в качестве выделенной (п +1) - ой траектории можно рассматривать любую из совокупности {0^}, г — 1,...,п + 1. Соответствующие различным сочетаниям произведения и суммы наиболее просто вычислять последовательным

~(' ' ~ 1 п+1 ~ (■)

через АЩ'. Пусть Ап — АШ' ■

=1

Теорема 2. Если выполняются условия леммы 1, причем р (Кр (а)) < 1 — е, (е > 0) У а, то ЕАп — Ап и ВА*п < БАп <

Доказательство. Доказательство строится путем прямого осреднения и оценок на основе сформулированных выше утверждений и условной независимости траекторий О1,..., 0п+1.

Оценку Ап естественно называть комбинированной. Отметим, что величина БАп может быть существенно меньше величины БАп щи Ба ^ 1 вследствие слабой коррелированное™ реализаций произведений в Ап для фиксированных значений в.

Аналогичная несмещенная элементарная оценка на основе леммы 1 для рп к Е \_.1'(а)/.1 (а)] строится по формуле

п (_1)я( о + 1) я

Ап — ?(0п+1,а) ■ ?(0п+2,а)^ У J,+2 П ,а) — -7о)

8=0 и0 к=1

с использованием п + 2 условно независимых траекторий. Более эффективная комбинированная оценка

п+1 (■)

Ап — Ап'/(п + 1) строится то аналогии с Ап для фиксированного значения £'(0п+2, а).

=1

3 Численные результаты

4.1. Для численного тестирования был рассмотрен односкоростной процесс переноса частиц с изотропным рассеянием (в том числе и после акта деления) в шаре радиуса Д — 7.72043 с параметрами: а — случайная величина, равномерно распределенная в интервале (0.7, 1.3);

— = 0.97; ^ = 0.03; у = 2.5; ас = 0; V— |у| — 1. а а

Как указано в [2], такой шар при а = 1с достаточной степенью точности критичен. Соответствующая преци-

А

приближении [6], которые дали значение \d(a = 1) = 0.000000 .... С использованием этого приближения путем моделирования значений а были получены следующие тестовые статистические оценки:

EAd = -0.00103, DAd = 0.00022.

Среднеквадратическая погрешность имеет здесь порядок последнего десятичного разряда. Для построения эффективного алгоритма метода Монте-Карло на основе теоремы 2 в сформулированную модель было введено поглощение с постоянным неслучайным коэффициентом acjv, который приводит к замене А I—у А — acjv Уа. Отметим, что такой прием является универсальным и может существенно повысить эффективность весового метода, исключая необходимость ветвления моделируемых траекторий, как это видно из дальнейшего.

Используя уравнение переноса (см., например [1]), можно сделать замену af — af + ас, v — vijfl(af + ac). В численном эксперименте моделировался процесс переноса с константами: а* = а\, а* = af + ас, v* = 1. Вспомогательные веса при этом определяются формулой Qn = {vaf/(af + ас))"*, где п\ — число делений, предшествующих данному столкновению [1]. Было нспользовано значение ас = 0.059, для которого vo f I (af + а с) < 1 и тем самым [1] р(Кр) < 1 Уст. Плотность распределения первых столкновений взята в виде

f(r,t)=2texp(-2t)fd(r), t> 0, r = |r| < R, (4.1)

где fd(r) = Сsin(^r)/r — диффузионное приближение к пространственной характеристической функции

для а = 1, ю = 0.3739866 [6]. Полагали также h(r, v) = sin(«(p)r)/г, где (см. [6])

«(,) = + Vaf )Р

R(as + vaf) + .

Расчеты показали, что использование таких функциональных параметров алгоритма существенно улучшает сходимость в (2.3) сравнительно с /¿(r) = (|-kR3) и даже с вариантом, в котором f (r,t)

определяется формулой (4.1), a h(r, v) = 1.

Нетрудно проверить, что для сформулированных выше характеристик вычислительной модели выполняются условия теорем 1, 2 и она тем самым позволяет построить оценку значений EA и DA. Для построения статистических оценок этих величин было реализовано 16 • 109 случайных значений оптической плотности а; моделирование траекторий осуществлялось параллельно на восьми процессорах Intel Core i7-3930K. При этом полагали п = 5, то есть для оценки EA использовади шесть, а для оценки EA2 — семь условно-независимых траекторий частиц. Результативные оценки определялись на основе анализа установления (с учетом вероятностной погрешности) двух старших значащих цифр при возрастании параметров tus. Таким образом получены следующие результативные оценки: EA = 0.00104 ± 0.00001 и DA = 0.00020 ± 0.00007 для t = 10 s = 4. В качестве погрешности здесь приведено среднеквадратическое отклонение. Эти оценки вполне соответствуют диффузионным, достаточная точность которых проверена лишь для а = 1.

Отметим, что полученные результаты переносятся на случай многоскоростного процесса переноса с анизотропным рассеянием соответственно [2].

4.2. Для проведения более сложных тестовых расчетов рассматривалась сферически симметричная среда с кусочно-постоянной случайной плотностью р = р(г) в шаре радиуса R = 7.72043 с макроскопическими

(0) (0) (0) сечениями рак ', pal , pa^ , где

а(0) = 1, а(0) = 0.97, af = 0.03, ^ = 2.5,

Для построения реализации среды шар делится па т одинаковых по объему сферических слоев, в каждом слое случайная величина р = pi выбирается независимо и равномерно на отрезке [1 — е, 1 + е].

Для такой системы в таблице 4.1 представлены результаты расчетов величины Ek из [7] с помощью специальной гомогенизации (то-есть путем осреднения реализаций р(г) с весом "диффузионной плотности" д(г)) и улучшенного диффузионного приближения; отметим, что в [7] рассматривалась величина к = р(К) — коэффициент размножения по поколениям столкновений. В таблице 4.1 приведены также оценки величин EA, DA с использованием этого приближения. Даны также основные здесь оценки, полученные методом Монте-Карло согласно изложенному в разделе 3. Представляет определенный интерес тот факт, что при т = 6 имеем EA > 0, хотя Ek < 1; результаты из [8] показывают, что при этом коэффициент размножения по поколениям делений keff > 1.

Таблица 4.1: Оценки EA и DA

Гомогенизация Метод Монте-Карло

т Ek EA VDA Ek EA VDA

1 0.9969 -0.00103 0.015 0, 9968 -0.00104 0.014

2 0.9981 -0.00063 0.012 0.9983 -0.00017 0.014

6 0.993 -0.00022 0.0069 0.9994 0.00020 0.007

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

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

[1] Марчук ГЛ.. Михайлов Г. А., Назаралиев М.А. и др. Метод Монте-Карло в атмосферной оптике. -Новосибирск: Наука, 1976, 283 с.

[2] Лотова Г. 3., Михайлов Г. А. Новые методы Монте-Карло для решения нестационарных задач теории переноса излучения // Журн. вычисл. матем. и матем. физики. 2002. Т. 42, №4. С. 569-579.

[3] Дэвисон Б. Теория переноса нейтронов. - М.: Атомиздат, 1960, 514 с.

[4] Марчук ГЛ. Методы расчета ядерных реакторов. М.: Госатомиздат, 1961.

[5] Larmier С., Zoia A., Malvagi F., Dumonteil Е., Mazzolo A. Neutron multiplication in random media: Reactivity and kinetics parameters // Annals of Nuclear Energy. 2018. V. 111. P.391-406.

[6] Романов Ю. А. Точные решения односкоростного кинетического уравнения и их использование для расчета диффузионных задач // Исследование критических параметров реакторных систем. М.: Атомиздат, 1960. С. 3-26.

[7] Лотова Г.З., Михайлов Г. А. Методы Монте-Карло для оценки вероятностных распределений параметров критичности процесса переноса частиц в случайно возмущенной среде // Журн. вычисл. матем. и матем. физики. 2018. Т. 58. № 11, С. 1900-1910.

[8] Лотова Г.З., Михайлов Г.А. Моменты параметров критичности процесса переноса частиц в случайной среде //Журн. вычисл. матем. и матем. физики. 2008. Т. 48. № 12, С. 2225-2236.

Михайлов Геннадий Алексеевич — чл. корр. РАН, Советник РАН Института вычислительной математики и математической геофизики СО РАН;

Новосибирский государственный университет;

e-mail: gam@sscc.ru;

Лотова Галия Зуфаровна — к.ф.-м.н., старший научный сотрудник Института вычислительной математики и математической геофизики СО РАН;

Новосибирский государственный университет;

e-mail: lot@osmf.sscc.ru;

Дата поступления — 1 июня 2019 г.