Математическое моделирование переноса излучения и переноса нейтронов с учетом процессов в сплошных средах
На правах рукописи
Аристова Елена Николаевна МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПЕРЕНОСА ИЗЛУЧЕНИЯ И ПЕРЕНОСА НЕЙТРОНОВ С УЧЕТОМ ПРОЦЕССОВ В СПЛОШНЫХ СРЕДАХ 05.13.18 – математическое моделирование, численные методы и комплексы программ
Автореферат диссертации на соискание ученой степени доктора физико-математических наук
Москва – 2009 год
Работа выполнена в Институте математического моделирования РАН
Официальные оппоненты:
член-корреспондент РАН, профессор Суржиков Сергей Тимофеевич доктор физико-математических наук, профессор Трощиев Виталий Ефимович доктор физико-математических наук, профессор Шагалиев Рашит Мирзагалиевич
Ведущая организация:
Физический институт им. П.Н.Лебедева РАН 1100 часов на заседании
Защита состоится « 9 » апреля 2009г. в диссертационного совета Д.002.058.01 при Институте математического моделирования РАН по адресу: Москва, 125047, Миусская пл, д. 4а.
С диссертацией можно ознакомиться в библиотеке ИММ РАН.
Автореферат разослан « » 2009г.
Ученый секретарь диссертационного совета Д.002.058.01, доктор физико-математических наук Змитренко Н.В.
ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ
Актуальность темы Решение уравнения переноса является одной из наиболее трудоемких частей моделирования задач, определяющую роль в которых играет перенос излучения и/или перенос нейтронов. Это связано с большой размерностью задачи, т.е. большим количеством переменных, от которых зависит функция распределения (фазовые переменные и время).
Построение эффективных методов решения уравнения переноса в многомерных геометриях является актуальной задачей при математическом моделировании динамических задач высокотемпературной радиационной газовой динамики (ВРГД) для различных приложений от взрыва сверхновой до инерциального термоядерного синтеза, в которых важнейшую роль в перераспределении энергии в системе имеет перенос собственного излучения плазмы.
Исследование саморегулируемых нейтронно-ядерных процессов в активных зонах быстрых реакторов приводит к необходимости динамического моделирования процессов переноса нейтронов совместно с выгоранием, реакторной кинетикой и управлением.
Трудности решения уравнения переноса помимо большой размерности связаны с несколькими факторами:
– построение аппроксимации дифференциального оператора в левой части уравнения связано с дилеммой точность-монотонность. По теореме Годунова любые линейные разностные схемы с порядком аппроксимации выше первого немонотонны, а во многих практически важных случаях еще и неположительны, что значительно снижает качество численного решения.
– весьма существенной проблемой является учет спектральной зависимости излучения (или частиц) из-за наличия резонансных областей коэффициентов поглощения, при этом коэффициенты в соседних точках спектра могут отличаться на несколько порядков.
– интеграл рассеяния в правой части уравнения приводит к итерационному процессу решения уравнения переноса, сходимость которого ухудшается, если индикатриса рассеяния содержит особенность преимущественного рассеяния вперед.
– во многих физических приложениях уравнение переноса нужно решать совместно с другими уравнениями, такими как уравнения газовой динамики в задачах переноса излучения, или уравнения выгорания и реакторной кинетики в задачах переноса нейтронов.
Объединенная система уравнений может обладать сильной нелинейностью, хотя само уравнение переноса линейно относительно своих переменных. Физически наглядно это для прохождения излучения через среду: поглощаясь, излучение меняет температуру среды, а измененная температура, в свою очередь, значительно изменяет коэффициенты поглощения.
Таким образом, взаимодействие излучения с веществом является нелокальным и нелинейным.
– при решении задач переноса нейтронов возможна постановка задачи на собственные значения, при решении которой в обычно используемых методах на итерационный процесс, связанный с рассеянием и делением, накладывается итерационный процесс нахождения собственного значения и собственной функции. Если при этом необходимо найти критическую сборку, т.е. определенное значение собственного числа, то возникает итерационный процесс третьего уровня, так что становится необходимым многократное решение уравнения переноса, что вычислительно очень дорого.
Создание эффективных численных методов для преодоления этих трудностей при моделировании переноса излучения и нейтронов с учетом процессов в сплошных средах очень актуально.
Цель и задачи исследования Настоящая работа посвящена разработке эффективных численных методов решения уравнения переноса в рамках квазидиффузионного подхода. Метод квазидиффузии (В.Я.
Гольдин. Квазидиффузионный метод решения кинетического уравнения // ЖВМ и МФ, 1964, т.4, № 6, с.1078-1087) заключается в постепенном понижении размерности используемых уравнений. На первом этапе происходит усреднение уравнения переноса по угловым переменным, в результате которого получается многогрупповая система уравнений квазидиффузии (КД). На втором – усреднение по спектру, что приводит к эффективной одногрупповой системе уравнений квазидиффузии. При внешнем усложнении подхода метод квазидиффузии позволяет решить некоторые из вышеперечисленных проблем. Во первых, метод квазидиффузии нелинеен, так как вводит дробно-линейные функционалы для вычисления компонент тензора квазидиффузии. Во-вторых, уравнения квазидиффузии выражают собой законы сохранения, поэтому консервативность получается автоматически.
В-третьих, при умеренной анизотропии рассеяния метод КД позволяет выразить главную часть рассеяния внутри группы через групповые скалярный и векторный потоки, что обеспечивает быструю сходимость итераций по рассеянию. В-четвертых, введение одногрупповой системы уравнений КД позволяет эффективно объединять эту систему с уравнениями, описывающими другие физические процессы внутри системы и учесть их взаимное влияние друг на друга. Для задачи на нахождение собственных значений и/или критических параметров системы значительно сокращается общее число итераций.
Разработка эффективных численных методов решения уравнений на всех этапах метода квазидиффузии в многомерных геометриях и приложения разработанных методов к математическому моделированию задач атмосферной радиации, управляемого термоядерного синтеза (УТС), а также к моделированию саморегулируемых режимов в активных зонах быстрых реакторов составляют содержание данной работы.
Методы исследования Методы работы основаны на построении разностных схем для дифференциальных уравнений в частных производных, построении методов решения полученных разностных уравнений и методов ускорения итераций. Проводится сопоставление численных решений с точными решениями, там, где они существуют, а также сравнение результатов математического моделирования с результатами других авторов и натурными экспериментами.
Научная новизна полученных в диссертации результатов состоит в следующем:
• Построены характеристическая и консервативно-характеристическая схемы решения уравнения переноса в собственных характеристических переменных для случая двумерной цилиндрической геометрии, учитывающие структуру логарифмических разрывов решения;
• Предложен аналог монотонной разностной схемы для аппроксимации несамосопряженных уравнений квазидиффузии, и построена комбинированная разностная схема решения уравнений квазидиффузии, сочетающая схему более высокого порядка точности в областях гладкости решения и аналог монотонной на контактных разрывах;
• Для несамосопряженной системы эллиптических разностных уравнений предложен нелинейный метод ускорения итераций решения эллиптических систем;
• Предложен метод учета сильной анизотропии рассеяния;
• Создан программный комплекс LATRANT для расчета задач газовой динамики при существенной роли собственного излучения плазмы;
• На основе разработанных математических моделей и созданных комплексов программ:
– решен ряд задач атмосферной радиации при наличии сильно анизотропного рассеяния с особенностью преимущественного рассеяния вперед на аэрозолях и в облаках;
– получены прецизионные расчеты теплового баланса атмосферы Земли для рэлеевского рассеяния на атмосферных газах при одновременном применении метода лебеговского усреднения по частотам, развитого А.В.Шильковым;
– проведено сравнение расчета задач УТС при учете излучения в многогрупповом приближении с трехтемпературной моделью плазмы, которое показало, что в трехтемпературной модели центральная область сжатия мишени имеет запаздывающую динамику по сравнению с многогрупповым приближением;
– проведено полномасштабное математическое моделирование экспериментов, проводимых в ТРИНИТИ, на лазерных установках PALS и LIL по взаимодействию мощных пучков лазерной энергии с пористыми средами;
– исследуются саморегулируемые нейтронно-ядерные режимы в активных зонах быстрых реакторов для создания реакторов нового поколения.
Теоретическая и практическая ценность результатов диссертации заключается: 1) в разработке эффективных численных методов решения уравнения переноса и уравнений квазидиффузии при наличии сильной анизотропии рассеяния совместно с уравнениями, отвечающими за другие физические процессы в полной системе;
2) в подробном исследовании задач УТС на основании предложенных моделей;
3) в предложении рекомендаций по формированию активных зон быстрых реакторов нового поколения, обладающих повышенными экономичностью и безопасностью по нейтронно ядерным процессам.
Основные публикации По теме диссертации опубликовано 54 работы. Основное содержание диссертации отражено в статьях [1-31].
Достоверность результатов диссертационной работы определяется их верификацией при разнообразном тестировании, включающем сравнение с точными решениями (при их наличии), сравнением с результатами экспериментов и расчетами по другим моделям, четким физическим смыслом полученных результатов и согласованностью с современными представлениями о предмете исследования.
Апробация результатов диссертации Результаты исследований, приведенных в диссертационной работе, были представлены и обсуждались на Всероссийских и Международных конференциях:
• Международный симпозиум «Численные методы решения уравнения переноса», Москва, май 1992;
• Fourth International Aerosol Symposium, St-Petersburg, July 1998;
• Международная конференция “Физика атмосферного аэрозоля”, Москва, апрель 1999;
• Joint International Conference on Mathematical Methods and Supercomputing for Nuclear Applications: Saratoga Springs, New York, October 1997;
Madrid, Spain, October 1999;
• Международная конференция “Математические идеи Л.П.Чебышева и их приложение к современным проблемам естествознания”, Обнинск, май 2002;
• XXVII European Conference on Laser Interaction with Matter ECLIM-2002, Москва, сентябрь 2002;
• VI International Congress on Mathematical Modeling, Нижний Новгород, сентябрь 2004;
• Международная конференция “Математика. Компьютер. Образование”, Дубна, январь 2004;
Дубна, январь 2006;
Дубна, январь 2008;
• XXXIII Международная (Звенигородская) конференция по физике плазмы и УТС, Звенигород, февраль 2006;
• International Congress on Advances in Nuclear Power Plants, Nice, France, May 2007;
• The 20th International Conference on Transport Theory, Obninsk, Russia, July 2007;
• 5th International Conference on Inertial Fusion Sciences and Applications, Kobe, Japan, September 2007.
Реализация и внедрение результатов работы Работа выполнялась в рамках научных планов Института математического моделирования РАН, проектов Российского фонда фундаментальных исследований, проектов МНТЦ, договоров с Физико-энергетическим институтом (Обнинск), Научно-исследовательским институтом атомных реакторов (Димитровград), РНЦ ”Курчатовский институт”.
Научные положения диссертации и разработанные на их основе методики, алгоритмы и программные комплексы использовались для совместных исследований в следующих организациях: Физико-энергетическом институте (Обнинск), Научно-исследовательском институте атомных реакторов (Димитровград), РНЦ “Курчатовский институт”, Физическом институте РАН.
Личный вклад соискателя В список положений, выносимых на защиту, включены результаты и выводы, в которых вклад соискателя был основным.
Структура и объем диссертации Диссертация состоит из введения, пяти глав, заключения, списка литературы, включающего 268 наименований, и изложена на 288 страницах.
ОСНОВНОЕ СОДЕРЖАНИЕ РАБОТЫ
Введение содержит обзор работ, относящихся к теме исследования и посвященных преодолению основных трудностей, возникающих при численном решении уравнения переноса;
обосновывается актуальность темы, выбор метода квазидиффузии как основного инструмента построения математической модели;
формулируется система уравнений ВРГД в рамках квазидиффузионного подхода, а также определяются вычислительные задачи, возникающие при решении этой системы уравнений при расщеплении по физическим процессам. Сформулированы цели проведения исследований, новизна диссертации и научные положения, выносимые на защиту. Приведены сведения об апробации работы и структуре диссертации.
Метод квазидиффузии заключается в постепенном понижении размерности задачи так, что на последнем этапе становится возможным эффективное объединение следствий уравнения переноса с уравнениями сплошной среды, описывающими другие физические процессы, например, с уравнениями газовой динамики для излучения или реакторной кинетики и выгорания для переноса нейтронов. Усреднением многогрупповых уравнений переноса по углам получается многогрупповая система уравнений квазидиффузии, замкнутая при помощи дробно-линейных функционалов, вычисляемых по решению уравнения переноса. Усреднение полученной многогрупповой системы уравнений квазидиффузии по энергии приводит к эффективной одногрупповой системе уравнений квазидиффузии относительно полных плотности и потока излучения, которая может быть объединена с другими уравнениями сплошной среды, что решает главную задачу – учета нелокального и нелинейного взаимодействия излучения и вещества. Внешнее усложнение процедуры решения уравнения переноса позволяет разрешить ряд трудностей. Во-первых, использование многогрупповой системы уравнений квазидиффузии автоматически означает консервативность схемы. Во-вторых, выделение главной части в члене рассеяния внутри группы позволяет перенести ее налево для уравнений квазидиффузии, так что влияние остаточных членов становится опосредованным через тензор квазидиффузии. Итерационный процесс решения уравнений с такой иерархической структурой быстро сходится при условии малости производных Фреше от введенных дробно-линейных функционалов по решению. Таким образом, решение уравнения переноса сводится к ряду вычислительных задач: 1) решение уравнения переноса при известной правой части, усреднение полученного решения по углам и вычисление тензора квазидиффузии, 2) решение многогрупповой системы уравнений квазидиффузии и усреднение этой системы по энергии в одногрупповую систему, 3) совместное решение полученной одногрупповой системы уравнений квазидиффузии с другими уравнениями сплошной среды (уравнениями энергии для задач радиационной газовой динамики или выгорания и реакторной кинетики). Диссертация посвящена методам решения поставленных вычислительных задач и моделированию физических задач на основе разработанных алгоритмов.
Глава 1 посвящена построению и сравнению двух схем решения уравнения переноса в собственных характеристических переменных. Метод итераций источника, обычно используемый при решении уравнения переноса при наличии рассеяния (и деления), сводится к тому, что правая часть вычисляется по решению, полученному на предыдущей итерации, поэтому при построении аппроксимации дифференциального оператора, стоящего в левой части уравнения, можно считать правую часть известной. Тогда уравнение переноса можно записать в виде:
+ = Q, (1) z + r + z r r = sin sin, [0, ], [0,2 ].
r = sin cos, z = cos, Переход к собственным характеристическим переменным позволяет записать уравнение, содержащее только две пространственные производные:
+ = Q, cos + sin (2) z s С формальной математической точки зрения связь r и с переменными s и h осуществляется по формулам (рис.1) r = s 2 + h s = r cos или (3) h = r sin = arctg(h / s ) + (1 sign( s)) / Рис. В отличие от классического метода Владимирова, в котором угловая сетка получается в результате пересечения множества параллельных касательных к окружностям постоянного радиуса пространственной сетки с множеством этих окружностей, в работе предложено независимое построение базовой угловой сетки для каждой из таких окружностей, например, равномерной по µ=cos. На практике это означает отказ от метода длинных характеристик и переход к методу коротких характеристик. С одной стороны, это сокращает вычислительные затраты и позволяет разрешить структуру логарифмических разрывов, характерных для слоистых задач в цилиндрической и сферической геометриях. С другой стороны, это усложняет алгоритм обхода ячеек и несколько понижает порядок сходимости метода из-за многократных интерполяций.
В диссертации предложены два алгоритма решения уравнения переноса в собственных характеристических переменных: характеристический и консервативно-характеристический.
В характеристическом методе значение в точке пересечения характеристики, выпущенной назад из узловой точки с неизвестным решением, с ребром расчетной ячейки восстанавливается со вторым порядком точности интерполяцией по значениям в трех освещенных узлах, после чего решение в четвертом узле находится интегрированием вдоль отрезка характеристики внутри ячейки. Консервативно-характеристический метод помимо значений в узлах использует интегралы от решения вдоль ребер ячейки (также называемые потоками), при этом квадратичная или псевдоквадратичная интерполяция на каждом ребре строится независимо от других ребер, что позволяет правильно учесть выпуклость функции.
В рамках предложенной аппроксимации на ребрах задача перераспределения выходящих потоков по неосвещенным граням решается точно.
В ходе решения задачи необходимо знать значение функции распределения в точках пересечения выпущенной назад характеристики с цилиндром предыдущего радиуса. Для этого предложен итерационный алгоритм построения монотонного сплайна смешанного порядка: второго, если это обеспечивает монотонность, и первого в противном случае.
Коэффициенты полученного сплайна используются в дальнейшем при вычислении компонентов тензора квазидиффузии интегрированием по азимутальному и полярному углам.
Для задач в сферической и цилиндрической геометрии характерно наличие логарифмических разрывов, связанных с принципиально разным поведением решения на двух геометрически близких характеристиках, если одна из них только касается области с другими параметрами, а вторая проходит по этой области. На логарифмическом разрыве одна из односторонних производных решения по азимутальному углу обращается в бесконечность. В предложенных методах в области контактных разрывов проводится близкая к касательной вторая характеристика, при этом при проведении интерполяции и при интегрировании учитывается поведение функции на логарифмическом разрыве между этими характеристиками, пропорциональное корню смещенного косинуса азимутального угла.
Проведено сравнение предложенных методов по порядку сходимости на задачах, имеющих точное решение. Показано, что использование консервативного варианта характеристической схемы значительно повышает порядок сходимости при незначительном удорожании расчета. Учет логарифмической структуры разрывов также улучшает качество предложенной схемы. Порядок сходимости консервативно-характеристического метода на владимировской сетке длинных характеристик – второй, а для метода коротких характеристик – чуть выше первого.
Результаты Главы 1 опубликованы в работах [1-2].
Глава 2 посвящена созданию и тестированию всех методов и программ, необходимых при численном решении системы квазидиффузионных уравнений и объединении их с уравнением энергии для вещества, что включает в себя: 1) построение разностных схем для многогрупповой системы уравнений квазидиффузии, 2) разработку методов решения полученных разностных уравнений эллиптического типа, 3) усреднение уравнений в эффективную одногрупповую систему, 4) решение объединенной системы усредненных уравнений квазидиффузии с уравнением энергии для вещества.
В расчетной (r-z) области вводится матрично упорядоченная сетка из четырехугольных ячеек. Построение разностной схемы для уравнений квазидиффузии (в заданной группе по энергии) ведется интегро-интерполяционным методом в рамках одной ячейки. Это упрощает ситуацию около контактных границ. В шаблон входят девять величин: значения плотности излучения и нормальные проекции потока в серединах четырех ребер, а в центре ячейки только плотность излучения (рис.2).
(i+1,j+1) Uj+ Wj+ Ui+ * W (i,j+1) * i+ U * (i+1,j) Ui * Wi *U j Wj (i,j) Рис.2. Разностная ячейка.
Ранее В.Я.Гольдиным и А.В.Колпаковым была предложена немонотонная аппроксимация уравнений квазидиффузии, основанная на кусочно-постоянной аппроксимации при вычислении контурных интегралов. Источников немонотонности у этой схемы два: во-первых, расширенный шаблон, так что полученная схема не удовлетворяет принципу максимума даже для самосопряженной задачи теплопроводности, во-вторых, несамосопряженность уравнений квазидиффузии в общем случае. Для преодоления первого источника немонотонности в работе предложена схема с уменьшенным шаблоном, для которой выполняется принцип максимума для самосопряженной задачи. Второе улучшение свойств монотонности достигается приведением тензора квазидиффузии к собственным осям в середине расчетной ячейки в плоскости r-z, что в силу предполагаемой непрерывности компонент уменьшает недиагональные компоненты тензора квазидиффузии на сторонах расчетной ячейки. Построение аналога монотонной схемы в повернутой системе координат приводит к значительному улучшению свойств монотонности схемы при полном учете несамосопряженности задачи.
Полученная система разностных уравнений решается методом µ-- прогонки, эквивалентной переносу граничных условий на стороны расчетной ячейки. При этом для нахождения µ и получается нелинейный быстросходящийся процесс, а для – линейный и довольно медленный. В диссертации предложен метод нелинейного ускорения этих итераций, основанный на введении отношений, вычисляемых в противоположных направлениях. Метод содержит итерационный параметр. Анализ сходимости метода при различных значениях итерационного параметра позволил предложить метод подстройки итерационного параметра под оптимальное значение в ходе итерационного процесса. Для несамосопряженной системы уравнений число необходимых для сходимости итераций пропорционально размерности разностной сетки по одному направлению.
Решение многогрупповой системы уравнений квазидиффузии усредняется в эффективную одногрупповую систему уравнений квазидиффузии для дальнейшего объединения полученной одногрупповой системы уравнений с уравнением энергии для вещества. При этом обменный член энергией между излучением и веществом выражается через произведение среднего коэффициента поглощения на суммарную по спектру плотность излучения и полную излучательную способность вещества. Эффективная линеаризация для применения метода Ньютона требует введения разностного аналога производной Фреше от усредненного коэффициента поглощения по температуре. Эта производная складывается из двух частей: первая отвечает за локальное просветление вещества при подъеме температуры, вторая, существенно нелокальная, отвечает за перестройку спектра на временном шаге вследствие изменения температур на расстояниях порядка длины пробега в окрестности данной точки. Эта производная вычисляется разностным образом на двух последовательных итерациях решения многогрупповой системы уравнений переноса с квазидиффузией на заданном временном шаге.
Предложенные методы тестировались на двух широко известных одномерных задачах Флека и их двумерных обобщениях. Первая задача Флека показала важность введения эффективной температуры усреднения при подготовке групповых констант для получения правильной скорости фронта волны изучения, вторая, обладающая контактной границей, показала важность введения производной Фреше от усредненного коэффициента поглощения (вплоть до численной разрешимости задачи).
Была рассмотрена задача, построенная на базе второй задачи Флека, о проникновении внешнего изотропного излучения с температурой 1 кэВ в трубу с внешним радиусом 1.1 см, внутренним – 1 см и длиной 5 см. Оптические характеристики внутренности трубы соответствуют умеренному коэффициенту поглощения в задаче Флека, а стенки обладают большим коэффициентом поглощения:
27 / T 3 (1 e ) при 0 r 1 см, ( ) = 4 (4) 10 / T 3 (1 e ) при 1 см r 1.1 см, (T ) = 8.11103 T. (5) Данная задача характерна тем, что фронт излучения в основной части трубы перпендикулярен контактной границе. На контактной границе происходит изменение направления движения волны излучения от распространения вдоль оси z для внутренности трубы к распространению по радиусу для стенок трубы. Это приводит к большим значениям недиагонального компонента тензора квазидиффузии около контактной границы (рис. 3 б), которыми нельзя пренебречь для построения аналога монотонной схемы. Использование аналога монотонной схемы с поворотом координат в плоскости r-z позволяет провести расчет этой задачи.
а) Drr б) Dzr в) Dzz Рис.3. Компоненты тензора квазидиффузии в восьмой группе, где находится максимум падающего планковского излучения. Справа от каждого рисунка приведена шкала линий уровня для момента времени 500 пс.
На рис. 4 приведены значения усредненного коэффициента поглощения по температуре для двух моментов времени, а также величины производных Фреше. Видно, как области с большим отрицательным значением этой производной (темные области), соответствующие просветлению вещества с ростом температуры, соседствуют с областями с большим положительным значением производной, что соответствует сдвигу спектра в низкоэнергетичную область.
a) t=0.25 нс б) t=0.5 нс в) t=0.25 нс г) t=0.5 нс Рис.4. Усредненный коэффициент поглощения (а,б) и производная Фреше по температуре от него (в, г) для моментов времени 250пс и 500 пс. Справа от каждого рисунка приведена шкала линий уровня.
Температурные поля для пяти моментов времени приведены на рис.5.
t=0.5 нс t=1 нс t=2. нс t=3 нс t=4 нс Рис.5. Динамика температуры.
Результаты Главы 2 опубликованы в работах [3-11].
Глава 3 посвящена методу учета сильной анизотропии рассеяния, которая характерна для задач переноса излучения в облаках и в присутствии аэрозолей в атмосфере, а также для ряда других задач от защиты реакторов до медицинских приложений использования лазера.
Ранее для умеренной анизотропии рассеяния было предложено в члене рассеяния выделять нулевой, первый и второй моменты разложения индикатрисы рассеяния по полиномам Лежандра и остаток. При использовании метода квазидиффузии главные члены рассеяния в уравнении переноса выражаются через плотность и поток излучения, а в уравнениях квазидиффузии соответствующие члены переносятся справа налево, обеспечивая значительное ускорение сходимости итераций по рассеянию, поскольку при таком преобразовании уравнения квазидиффузии зависят от рассеяния опосредовано через коэффициенты квазидиффузии. Однако наличие в индикатрисе рассеяния сильной особенности преимущественного рассеяния вперед делает такое преобразование недостаточным: во-первых, при этом очень плохо затухают коэффициенты разложения индикатрисы по полиномам Лежандра, и остаток разложения велик, и, во-вторых, в этом случае восстановление индикатрисы по первым трем членам не гарантирует положительности основной части члена рассеяния. В диссертации предложено в этом случае выделять -образную особенность индикатрисы рассеяния с малым носителем в окрестности углов прямого пролета без рассеяния. Оставшееся гладкое продолжение индикатрисы частично разлагается по полиномам Лежандра, а остатки интеграла рассеяния после выделения главных частей регулярной и сингулярной компонент индикатрисы преобразуются в дробно-линейные функционалы. Для вычисления интегралов в этих функционалах предложена замена переменных, позволяющая учесть другие возможные особенности индикатрисы рассеяния. В предложенном методе аналогичным образом предложено учитывать анизотропию отражения в граничных условиях.
Сходимость метода тестировалась на задаче рассеяния в однородном плоском слое 0z10, граничащем с вакуумом с двух сторон, с изотропным источником q=1, полным сечением =1 и сечением рассеяния s=с. Для индикатрисы рассеяния Хеньи-Гринстейна 1 g2 1 1+ g wHG ( µ ) = w(1) =, 2 (1 + g 2 2 g µ )3 / 2 2 (1 g ) была исследована зависимость числа итераций по рассеянию от параметров жесткости задачи g и c. Здесь g – параметр, отвечающий за особенность преимущественного рассеяния вперед. Чем ближе величина g к единице, тем ближе особенность к -особенности при µ1.
В обычных методах итераций источника наихудшая сходимость имеет место при отсутствии поглощения. Количество итераций по рассеянию в предложенном методе (последние четыре колонки) в сравнении с рядом других методов из работы А.В.Волощенко представлено в Табл. 1.
Наибольшее число итераций по рассеянию возникает для задачи средней жесткости (параметр g=0.99) при отсутствии выделения -особенности. В остальных случаях число итераций невелико и колеблется от 9 до 15. Ужесточение критерия сходимости ведет к двукратному росту числа итераций при замене на 2.
Таблица 1.
ik ik = 104 при изменении параметров g и c.
Число итераций при критерии сходимости max k i i Число итераций µ P1SAa N/Ab P1SA-c g c 0.95 0.99 0.999 1.
0.9 0.9 20-15 66-45 9 11 0.99 66-23 392-89 12 17 0.99 0.9 28 48 12 0.99 75 120 13-21 15 0.999 0.99 5-7 Далее в диссертации рассмотрен ряд задач атмосферной радиации. Прохождение излучения в атмосфере связано с двумя классами задач. Во-первых, прямые задачи, когда по известным распределениям веществ с заданными свойствами требуется найти поле излучения атмосферы или его (излучения) какие-либо интегральные характеристики. Во вторых, обратные задачи, т.е. задачи диагностики состояния атмосферы по прохождению сигнала (роль которого может играть и естественный источник – Солнце). В обоих случаях немаловажную роль играет сильно анизотропное рассеяние на аэрозолях и в облаках.
В диссертации предлагается метод корректного учета сильной анизотропии рассеяния, примененный для исследования рассеяния солнечного излучения в атмосфере, содержащей аэрозоли и облачные слои. Рассеяние в облаках и на аэрозолях характеризуется дельтаобразной особенностью преимущественного рассеяния вперед, т.е. большой долей рассеяния на малые углы (рис.6а). Рассеяние на малые углы существенно при наличии сосредоточенных источников, роль которых в задачах атмосферной радиации играет солнечное излучение.
Используется приближение плоской атмосферы, широко распространенное для климатических расчетов. Это позволяет не учитывать азимутальную зависимость интенсивности падающего и рассеянного света I, сведя ее к равномерному распределению относительно азимутального угла (), что уменьшает число независимых переменных до одной пространственной переменной () и одной угловой (µ) – косинусу угла с нормалью, направление которой выбрано от верхней границы атмосферы к земле. Однако разность азимутальных углов между направлением распространения и рассеяния света входит как (1 µ )(1 µ ) cos, от которого зависит 2 параметр в косинус угла рассеяния = µµ + индикатриса рассеяния w, и по этому углу ведется интегрирование в интеграле рассеяния.
Переход к многомерным задачам в рамках данной методики сведется к увеличению независимых переменных в уравнении переноса и связанных с этим проблемам, не ограничивая применимость данного метода учета анизотропного рассеяния.
-особенности Рис.6а. Индикатриса рассеяния для модели водяного Рис.6б. Выделение в облака с широким распределением частиц по индикатрисе рассеяния. Пунктиром размерам, обладающая сильной особеннос- нарисована регулярная компонента тью преимущественного рассеяния вперед.. wreg().
Стационарное уравнение переноса в приближении плоской атмосферы для света заданной частоты имеет вид:
I (, µ ) 1 + ( + ) I (, µ ) = w (, )I (, µ ) d d µ, µ (6) 1 где – коэффициент поглощения, а – коэффициент рассеяния.
В предлагаемом методе интенсивность излучения I представляется в виде суммы трех компонент:
1) I0 – нерассеянное солнечное излучение, в том числе зеркально отраженное от поверхности земли, 2) I1 – излучение, рассеянное на малые углы, в том числе с многократным рассеянием, 3) I2 – все оставшееся излучение.
На верхней границе атмосферы задается падающее солнечное излучение в виде функции угла: I=q·(µ-µ0). На нижней границе атмосферы, т.е. на поверхности земли с альбедо r, задается условие отражения, при этом предполагается, что доля излучения отражается зеркально, а доля 1- – диффузно. Метод позволяет учитывать и более сложные варианты граничных условий отражения, но за неимением соответствующих данных пришлось ограничиться комбинацией зеркального и диффузного отражения.
Решение задачи для компоненты I0 выписывается аналитически в виде затухающей экспоненты от текущей оптической толщины, деленной на косинус угла падения.
Учет однократного рассеяния компоненты I0 в компоненте I1 проводится полуаналитически с высокой степенью точности, для многократного рассеяния этой компоненты решения использованы алгоритмы ускорения итераций на основе потоковой квазидиффузии, которая является естественной здесь, поскольку решение для этой компоненты распадается на две ветви с существенными значениями в окрестности углов падения и зеркального отражения ±µ0.
Уравнение переноса для компоненты I2 содержит решения для I0 и I1. Для его решения используется метод квазидиффузии, модифицированный для учета сильной анизотропии рассеяния.
Отметим характерные черты предложенного метода.
1. Вычисление интеграла рассеяния.
a) Метод позволяет использовать максимально полную информацию об индикатрисе рассеяния при ее наличии, не ограничиваясь знанием только нескольких первых коэффициентов разложения индикатрисы рассеяния по полиномам Лежандра.
b) Для корректного учета рассеяния на малые углы индикатриса рассеяния делится на гладкую регулярную часть с областью определения [-1,1], и сингулярную, с малым носителем в окрестности µ=1: w()=wreg()+wsing() (рис.6б). Главная часть сингулярной компоненты индикатрисы рассеяния может быть представлена в виде -функции, подстановка которой в уравнение переноса приводит к эффективному уменьшению коэффициента рассеяния. При этом в правой части остается поправочный член рассеяния на малые углы, который можно не учитывать в тех областях фазового пространства, где рассеяние на малые углы несущественно (вдали от сосредоточенного источника).
c) Вычисление нулевого, первого и второго коэффициентов разложения регулярной компоненты индикатрисы по полиномам Лежандра используется для выделения главной части в интеграле регулярного рассеяния, что существенно ускоряет сходимость итераций по оператору рассеяния.
d) После выделения главных частей в регулярной и сингулярной составляющих индикатрисы рассеяния, соответствующие интегралы от остатков преобразуются в дробно-линейные функционалы, что также ускоряет сходимость итераций.
e) Вычисление интегралов в этих дробно-линейных функционалах может производиться либо напрямую, либо с заменой переменных µ/,,. Эта замена позволяет хорошо учитывать особенности индикатрисы рассеяния (wsing), но при этом интерполируется интенсивность излучения. В зависимости от поведения функций, имеется возможность интерполировать более гладкую из них, выбирая либо непосредственное вычисление интеграла, либо с заменой переменных.
f) Если выбран метод вычисления интеграла рассеяния от остатков индикатрисы с использованием замены переменных, сетка интегрирования по продуцируется сеткой, на которой задана индикатриса рассеяния (см. пп. а)), в силу однозначности связи и при фиксированных µ,µ/,[0,].
g) Член рассеяния на малые углы после выделения главного члена wsing в предположении дифференцируемости интенсивности излучения может быть приведен к Фоккер Планковскому виду, т.е. к диффузии в фазовом пространстве. Однако мы не используем это представление, поскольку с нашей точки зрения предположение о дифференцируемости функции распределения является очень сильным предположением, которое вряд ли справедливо в окрестности µ=0 и µ=±µ0.
2. Для ускорения сходимости итераций по интегралу рассеяния использован метод квазидиффузии, в котором наряду с уравнением переноса используются его макроскопические следствия – уравнения для плотности и потока излучения, замкнутые при помощи введения дробно-линейных функционалов. При этом интеграл рассеяния также представляется в виде, максимально использующем информацию о плотности и потоке излучения (см. пункт 1.с)). Скорость сходимости итераций по рассеянию такой расширенной системы определяется величиной производных Фреше введенных дробно линейных функционалов, и в рассмотренных задачах скорость сходимости оказывается весьма высокой.
3. Численная схема как для уравнения переноса, так и для уравнений квазидиффузии выписывается в предположении постоянства коэффициентов поглощения и рассеяния на каждом расчетном интервале. Нерассеянное солнечное излучение входит в виде экспоненциального источника в уравнения для компонент I1, I2, кроме того, в уравнения для I2 входит решение для компоненты I1. Численная схема учитывает экпоненциальную зависимость I0, а зависимость всех остальных членов по пространству предполагается кусочно-линейной.
Методика расчета поглощения и рассеяния солнечного излучения при заданной частоте была опробована на нескольких модельных задачах атмосферной радиации для вертикального и наклонного углов падения солнечного излучения с индикатрисой рассеяния, характерной для рассеяния в облаке с широким распределением частиц по размерам (рис. 6) и обладающей очень сильной особенностью преимущественного рассеяния вперед. Было выявлено влияние каждой из компонент решения на уширение солнечного пучка, а также отмечено возникновение разрывов решения на границах облака.
Была решена задача о рассеянии солнечного излучения для реальной модели атмосферы, состоящей из воздуха с рэлеевским законом рассеяния и ряда модельных аэрозолей (фонового стратосферного, континентального, городского, морского – рис. 10), при отсутствии поглощения. Рассеяние на аэрозолях также обладает особенностью преимущественного рассеяния вперед, соответствующие индикатрисы рассеяния были рассчитаны в книге (М.Я.Маров, В.П.Шари, Л.Д.Ломакина. Оптические характеристики модельных аэрозолей атмосферы земли. ИПМ им. М.В.Келдыша АН СССР, М.: 1989, 230с) и представлены на рис.7. Из предложенных для сравнения вариантов мы выбрали расчет для длины волны 0.55µ при нормальном падении солнечного излучения. Распределение аэрозолей соответствует профилю II рабочей группы «Стандартная радиационная атмосфера» (Перенос радиации в рассеивающих и поглощающих атмосферах. Стандартные методы расчета. Под ред. Ж.Ленобль. – Л.: Гидрометеоиздат, 1990, 230с.) и представлено в Табл. 2. Все остальные данные также соответствуют этому варианту. Расчет велся в физических переменных в слое 104км0км.
Таблица 2.
Распределение аэрозолей и оптические толщины для реальной модели атмосферы.
–– Стратосферный аэрозоль Городской |\ –– –––––––– =0.003 ––––––––– аэрозоль |\ =0.5 км- –– Континентальный |\ з =2.18·10-4 км- –– =1. |\ е –– аэрозоль |\ м Континентальный –– |\ л аэрозоль =0.0025 км-1 =0.1 км- –– |\ я –– |\ =0. =3.32·10-5 км- –– =0.025 |\ Морской –– |\ аэрозоль =0.025 км- –– |\ –– |\ =0. –– |\ ––––––––––––– + Рэлеевское рассеяние атмосферных газов –––––––––––– –– |\ =0. -104 км -30 км -20 км -12 км -2 км 0 () Для чисто рассеивающей атмосферы сохранение полного радиационного потока по высоте является точным следствием уравнения переноса. Данные о потоке при различных углах падения солнечного излучения при наибольшей оптической толщине атмосферы по рассеянию в сравнении с лучшими результатами рабочей группы «Стандартная радиационная атмосфера» приведено в Табл. 3. Предложенный метод сохраняет поток по высоте и при наклонном падении, что эффективно также увеличивает толщину атмосферы по рассеянию.
Задача IV. Индикатрисы рассеяния континентального, морского, городского и Рис.7.
стратосферного аэрозолей с увеличенным масштабом углов преимущественного рассеяния вперед.
Таблица 3.
Полный поток на длине волны 0.55 мкм для чисто рассеивающей атмосферы, содержащей атмосферные газы с рэлеевским рассеянием и ряд аэрозолей с сильно анизотропным рассеянием. Сравнение с результатами рабочей группы «Стандартная радиационная атмосфера» в случае, когда в нижнем двухкилометровом слое имеется городской аэрозоль с оптической толщиной по рассеянию, равной 1.
Косинус угла Результат предложенного Результаты рабочей группы падения метода солнечного Поток постоянен по z=30км z=2км z=0км излучения высоте µ0=1. 2.575 2.863 2.852 1. µ0=0.75 1.759 2.070 2.059 1. µ0=0.5 0.997 1.292 1.282 0. А.В.Шильковым с соавторами была разработана система кодов атмосферной радиации ATRAD, включающая восстановление коэффициентов поглощения атмосферными газами по параметрам линий из базы данных HITRAN и подготовку коэффициентов поглощения и рассеяния, усредненных по Лебегу. Вид уравнения переноса сохраняется при использовании лебеговских групп. Точность метода лебеговского усреднения на 400 эффективных группах совпадает с точностью методов типа ‘line-by-line’ при сокращении на четыре порядка количества используемых энергетических точек. Разработанная методика учета анизотропии рассеяния была включена в систему кодов атмосферной радиации ATRAD.
Результаты расчетов по предложенной методике при учете как поглощения всеми атмосферными газами, так и молекулярного рэлеевского рассеяния, представлены в Таблице 4.
Таблица 4.
Полное поглощение всеми атмосферными газами в линиях с учетом и без учета молекулярного рассеяния (Вт/м2).
Альбедо поверхности 0.2 Альбедо поверхности 0. =30 =75 =30 = Поглощение 178.2 73.8 200.4 81. Поглощение 175.6 71.0 197.5 79. и рассеяние Изменение потока по высоте для случая поглощения всеми атмосферными газами без учета и с учетом рассеяния приведено на Рис. 8. Результатов ‘line-by-line’ расчетов для рассеивающих атмосфер крайне мало.
На рис. 9 приведена зависимость скорости радиационного выхолаживания от давления, полученная методами ‘line-by-line’ разными авторами для чисто поглощающей атмосферы и по расчетам в системе ATRAD.
Таким образом, можно заключить, что метод лебеговского усреднения дает точность расчета методов ‘line-by-line’ при сокращении объема вычислений в 10000 раз.
Предложенный метод учета анизотропии рассеяния хорошо работает как в лебеговских обобщенных, так и в физических переменных. При объединении метода лебеговского усреднения с методом учета сильной анизотропии рассеяния удается получить прецизионные результаты расчетов по тепловому балансу атмосферы Земли с учетом рассеяния, которые практически недоступны для других методов расчета.
Результаты Главы 3 опубликованы в работах [12-17].
Рис.8. Поток солнечного излучения в стандартной летней атмосфере средних широт при молекулярном поглощении в линиях всеми атмосферными газами (сплошная линия) и при учете рассеяния (пунктир).
Рис.9. Скорость радиационного выхолаживания для стандартной летней атмосферы средних широт.
Глава 4 посвящена описанию программного комплекса LATRANT, созданного для математического моделирования задач инерциального термоядерного синтеза.
Программный комплекс учитывает перенос излучения в многогрупповом приближении и движение газа по улучшенной лагранжевой методике в двухтемпературном приближении в переменных r-z. Программный комплекс был создан на базе газодинамического комплекса ATLANT, созданного А.Б.Искаковым, и программного комплекса LATRA совместного решения уравнений переноса излучения и энергии вещества, созданного автором диссертации.
Уравнения газовой динамики в двухтемпературном односкоростном приближении для электронного и ионного компонент плазмы с учетом переноса излучения в лагранжевой системе координат имеют вид d + div u = 0, (7) dt W du + grad( pe + pi + ) = d, (8) a dt c d e + div We + ( pe + ) div u = Qie + Qr + Qe + Qz, (9) dt di + div Wi + ( pi + (1 ) ) div u = Qie + Qi. (10) dt Здесь:
Ti Te обменный член энергией между электронной и ионной Qie = компонентами плазмы;
I ( Z ) отвечает за кинетику ионизации, в этой формуле Z – Z F ( Z,, Te ) Qz = Z заряд иона, I(Z) – энергия ионизации, F ( Z,, Te ) – скорость ионизации;
член обмена энергией между излучением и веществом, Qr = ( U U )d Pl Pl U – спектральная плотность излучения, планковская равновесная плотность излучения 8 h U = 4 B = Pl exp(h kTe ) c3 (h – постоянная Планка, k – постоянная Больцмана, – постоянная Стефана-Больцмана), – спектральный коэффициент поглощения с поправкой на вынужденное переизлучение.
a Остальные обозначения универсальны. В дальнейшем мы пренебрегаем членом давления излучения в уравнении движения (2), поскольку оно существенно только для сверхвысоких температур.
Система (1) – (7) замыкается уравнениями состояния:
pe = pe (, Te ), pi = pi (, Ti ), e = e (, Te ), i = i (, Ti ).
На границах расчетной области задается либо условие непротекания, либо значения внешнего давления, а также значения тепловых потоков для электронного и ионного компонента.
Для спектрального описания задачи используется многогрупповое приближение, уравнение переноса в квазистационарном приближении в r-z геометрии имеет вид:
I p I p p 1 I p I p w + (U + sp ) I p = s ( ) I p d + Pl B p + z + r + p p p (11) r k c t z r z = cos, r = sin cos, = sin sin, d = sin d d, где p – индекс группы, wkp – индикатриса рассеяния в данной группе, sp – коэффициент рассеяния в группе, а по спектральному коэффициенту поглощения вычисляются три a средних групповых коэффициента поглощения по формулам:
p +1 p + Pl (T, ) = (T,, )U Pl (T, )d p U Pl (T, )d, (12) p p p +1 p + p U (T,, ) = (T,, )U Pl (, )d U Pl (, )d, (13) p p p + U Pl (, ) U Pl (, ) p + 1 = d d, (14) p R (T,, ) p (T,, ) p последний из коэффициентов используется в многогрупповых уравнениях квазидиффузии.
Здесь – эффективная температура усреднения. Используемая в расчетах база данных p DESOPLA (ФИ РАН) включает в себя таблицы групповых коэффициентов поглощения Pl, p p U, R для значений плотностей 10-5, 10-4, 10-3, 10-2, 0.1, 2.5, 20, и 100 г/см3, температур 2.6·10-5, 10-3, 10-2, 5.·10-2, 0.1, 0.5, 1, 3, 10 кэВ, и значений эффективной температуры :
0.5 кэВ, 1 кэВ, 2 кэВ, 3 кэВ для наиболее часто используемых в экспериментах по УТС веществ. Под первоначально понималось максимальное значение температуры в расчете, которое должно было либо априорно оцениваться, либо определяться из предварительного расчета. В дальнейшем была сделана модификация использования всей базы данных по, используя оценку величины по положению максимума пришедшего излучения.
Применение метода квазидиффузии для решения уравнения переноса позволяет эффективно учитывать нелинейное и нелокальное взаимодействие излучения с веществом.
Расщепление по физическим процессам приводит к следующей схеме расчета уравнений на одном временном шаге:
0. Запоминание всех необходимых нестационарных величин с предыдущего временного шага;
1. Расчет уравнений газовой динамики с искусственной вязкостью – разлет ячеек с дробными шагами по времени, величину которых диктует явный алгоритм газовой динамики;
вычисление вклада искусственной вязкости в уравнения энергии;
2. Расчет вклада электронной и ионной теплопроводности в уравнения энергии;
3. Решение групповых уравнений переноса, их усреднение по направлениям полета фотонов в квазидиффузионные многогрупповые уравнения;
4. Суммирование и усреднение групповых уравнений квазидиффузии по энергии в эффективную одногрупповую систему уравнений квазидиффузии;
5. Совместное решение двух уравнений энергии для электронного и ионного компонента плазмы и эффективной одногрупповой системы уравнений квазидиффузии;
6. Расчет энергетического баланса.
Для решения многогруппового уравнения переноса излучения совместно с уравнениями энергии для электронного и ионного компонент вещества используются методики, описанные в Главах 1 и 2.
Самым дорогим в вычислительном отношении является решение уравнения переноса.
Поэтому для удешевления расчета могут использоваться некоторые приближения, например, многогрупповое диффузионное приближение (не учитывающее анизотропии распространения излучения по углам) или еще более простая трехтемпературная модель, в которой излучение описывается планковской функцией с радиационной температурой.
На основании предложенного программного комплекса был решен ряд задач инерциального термоядерного синтеза. Как известно, гидродинамические неустойчивости и перемешивание препятствуют достижению давлений и температур, необходимых для начала термоядерной реакции. Поэтому физически очень важной задачей является достижение максимальной симметрии обжатия мишени. Увеличение количества лазерных пучков технически значительно усложняет каждый выстрел из-за проблемы синхронизации пучков и не позволяет разрешить проблему неоднородности из-за интерференционных явлений внутри самих высококогерентных пучков, их перекрытия и еще ряда технических и физических трудностей. Для улучшения симметрии сжатия мишени было предложено несколько подходов. Первый из них – это переход к мишеням типа «лазерный парник», в которых лазерное излучение поглощается и переизлучается в рентгеновском диапазоне стенками камеры, и уже это рентгеновское излучение сжимает мишень. При этом энергия лазера, идущая собственно на сжатие мишени, значительно уменьшается. Второй физической идеей было использование симметризующего предимпульса для создания плазменной короны до прихода основного лазерного импульса. Однако при этом часть энергии предимпульса переизлучается в мягком рентгене, что также ухудшает условия сжатия центральной области основным импульсом. Третий вариант увеличения симметрии обжатия мишени заключается в покрытии мишени малоплотной пеной. Одной из главных особенностей малоплотных сред является их способность сжиматься под действием ударной волны до плотностей в несколько раз более высоких, чем те же вещества с полной плотностью. Однако удаление от центра сферы зоны поглощения лазерного излучения приводит к снижению эффективности мишени. Перспективность использования слоя абсорбера-аблятора при сознательном уменьшении эффективности мишени, но при достижении устойчивости сжатия, вызвало волну как экспериментальных, так и теоретических исследований на основе различных математических моделей. Основными механизмами сглаживания неоднородностей являются диффузионная и радиационная теплопроводность. Для увеличения радиационной теплопроводности было предложено использовать либо покрытие пен тонкой пленкой из тяжелых атомов (золота или меди), или добавлять такие тяжелые атомы непосредственно в пену.
На основе созданного комплекса LATRANT:
1. Было проведено математическое моделирование мишеней типа «лазерный парник» изотропным внешним излучением с температурой 150эВ. Моделирование сферически симметричных задач в двумерной постановке показало, что все схемы хорошо поддерживают симметричность при сжатии и разлете мишени. Было проведено сравнение результатов сжатия и разлета модельной четырехслойной двухоболочечной мишени из стекла и дейтерий-тритиевой смеси между и внутри стеклянных оболочек при учете переноса излучения в многогрупповом диффузионном приближении и в трехтемпературном приближении. Было показано, что использование трехтемпературной модели дает запаздывающую динамику сжатия центральной области горючего по сравнению с многогрупповым расчетом (рис. 10) за счет преднагрева центральной области горючего высокоэнергетичными фотонами и большую температуру центральной области горючего.
Рис.10. Динамика сжатия центральной области DT горючего (слева) и электронная температура центральной области (справа).
2. Была рассчитана существенно двумерная задача о сжатии той же сферической мишени внешним изотропным излучением с температурой, меняющейся в зависимости от полярного угла. Результаты расчетов показывают, что наблюдается симметризация сжатия внутренней оболочки по сравнению с несимметричным сжатием внешней оболочки.
3. Было проведено исследование взаимодействия лазерного излучения с фольгами для дальнейшего сравнения с результатами покрытия фольг малоплотными пенами.
Плазменная корона фольги в три раза увеличивает радиус светящейся с тыльной стороны области по сравнению с радиусом фокального лазерного пятна. При этом выходящее излучение становится значительно более изотропным по радиусу.
4. Были рассчитаны две задачи о взаимодействии лазерного излучения с двумя мишенями, состоящими из пористого полиэтилена на алюминиевой подложке. В первой задаче эта мишень со стороны падения лазерного излучения была покрыта тонкой пленкой меди, а во втором случае кластеры меди были распределены равномерно в пене. Количество атомов меди на единицу площади во втором случае было в полтора раза меньше, чем в первом. Эффективность конверсии лазерного излучения в рентгеновское показана на рис.11.
Сравнение компонент баланса энергии показывает, что во втором случае до 80% лазерной энергии высвечивается. Эта задача показывает, что хотя доля собственно энергии излучения в этих задачах ничтожно мала, именно перенос излучения осуществляет обмен энергией между отдельными частями системы.
Рис.11. Распределение энергии по составляющим для задачи I (слева) и II (справа).
На рис.12 представлен поток излучения со стороны подложки этих двух мишеней.
Видно, что в задаче с кластерами меди существует момент времени (около 2 нс), после которого поток излучения нарастает за очень короткое время. Этот эффект регистрируется в экспериментах: для пены данной толщины (порядка 500 мкм) время появления интенсивного свечения в эксперименте составляло также 2 нс.
Рис. 12. Зависимость от времени интегрального потока излучения с обратной стороны мишени для первой (слева) и второй (справа) задач.
5. На установке PALS (Prague Asterix Laser System) была проведена серия экспериментов по взаимодействию излучения йодного лазера с плоскими пористыми мишенями.
Энергия лазерного импульса в эксперименте оценивалась в 150-175 Дж. Длительность импульса около 0.8 нс с шириной импульса на полувысоте 320 пс. Такие значения 5 1014 Вт/см2 на параметров соответствуют потоку лазерной энергии порядка поверхности мишени.
Использовались пористые мишени из триацетата целлюлозы на алюминиевой подложке толщины 5 мкм. Толщина пены около 400 мкм. При этом предполагались две основные плотности пористого вещества – 9.1·10-3 г/см3 и 4.5·10-3 г/см3:
=2.7 г/см -0.0405z-0.0400 см Al (толщина слоя 5мкм) 1=9.1·10-3 г/см3;
(толщина слоя 400мкм) -0.0400z 0.0000 см TAC 2=4.5·10-3 г/см Для ряда выстрелов использовалась пена с добавками 9.9% по массе атомов меди 1=9.1·10-3 г/см3;
TAC +9.9% Cu Первоначально для сравнения расчета с экспериментом в качестве поглощенной лазерной энергии была заложена энергия 130Дж. На рис.13 приведены распределения температур для моментов времени, соответствующих максимуму падающего излучения, и ближе к окончанию лазерного импульса.
Рис. 13. Электронная температура для пены ТАС (9.1 мг/см3, 3 гармоника лазера) на моменты времени 400, 500 и 600пс, и ионная температура на момент 600пс.
В эксперименте регистрировались показания рентгеновской электронно-оптической камеры (РЭОК) прихода гамма-квантов с энергией больше 1.5 кэВ. Фотографии РЭОК по горизонтальной оси дают временную развертку регистрируемого рентгеновского излучения с полным временем 2 нс по ширине фотографии, а по вертикальной оси – развертку по толщине мишени с полной толщиной 2000мкм по высоте фотографии. По этим данным вычислялись две скорости: скорость распространения рентгеновского фронта как касательная к нижней части цветного изображения в носике на границе с черным фоном и скорость гидротепловой волны по сечению этого рисунка в максимуме регистрируемого излучения. Пример фотографий с РЭОК приведен на рис.14.
a) выстрел № 28205: Elas=170,4 Дж;
2нс. b) выстрел № 28232: Elas = 163 Дж;
2нс.
Рис.14. Показания РЭОК для пены TAC (9.1 мг/см3, 3 гармоника лазера, толщина фольги 5 мкм) в экспериментах (а, б) на установке PALS.
Программный комплекс позволяет моделировать показания РЭОК. Результат такого моделирования приведен на рис.15. Скорости распространения рентгеновского фронта и гидротепловой волны, измеренные по приведенным данным, оказались значительно выше, чем соответствующие скорости, измеренные по экспериментальным данным. Заметим, что сравнение рис. 14 а) и рис. 14 b) показывает, что в экспериментах не достигается полной повторяемости результатов: если на левой картинке видно свечение алюминиевой подложки ко времени примерно 1.5нс, то на правой картинке этого свечения не наблюдается даже и к 2нс для почти одинаковой энергии выстрелов. В ФИ РАН по теории сильного взрыва были проведены оценки вложенной лазерной энергии, которые показали, что регистрируемые в эксперименте скорости ударных волн соответствуют значительно меньшим вложенным энергиям. Решение одномерных уравнений Максвелла, полученное численно также в ФИ РАН, показало, что на слоистой структуре пленок из ТАС отражается до 70% лазерной энергии. Поэтому было проведено исследование задач взаимодействия лазерного излучения с пенными структурами с меньшими значениями вложенной лазерной энергии. Результаты моделирования для пены 9.1 мг/см3 при облучении на третьей гармонике йодного лазера при различных энергиях приведено на рис.15. Аналогичные данные, полученные для пены 4. мг/см3 при облучении мишени на первой гармонике лазера в эксперименте представлены на рис.16, а при математическом моделировании – на рис. 17.
а) Eabs= 130 Дж;
0,6нс b) Eabs= 130 Дж c) Eabs = 50 Дж;
1нс d) Eabs = 50 Дж e) Eabs = 6 Дж, 1нс f) Eabs = 6 Дж Рис. 15. Результаты моделирования показаний РЭОК (слева) и энергетический баланс системы (справа) для пены TAC (9.1 мг/см3, 3 гармоника лазера, толщина фольги 5 мкм) в расчетах при разных энергиях выстрелов.
a) выстрел № 28256: Elas= 174 Дж ;
2нс b) выстрел № 28270: Elas = 177 Дж ;
2нс Рис. 16. Показания РЭОК для пены TAC TAC (4.5 мг/см3, 1 гармоника лазера, толщина фольги 5мкм) в экспериментах.
a) Eabs = 50 Дж ;
1нс b) Eabs = 50 Дж c) Eabs = 10 Дж, 1нс d) Eabs = 10 Дж Рис. 17. Результаты моделирование показаний РЭОК (слева) и энергетический баланс системы (справа) для пены TAC (4.5 мг/см3, 1 гармоника лазера, толщина фольги 5мкм) в расчетах при разных энергиях выстрелов.
Добавление тяжелых кластеров меди (9.9% по массе) в пену ТАС значительно усиливает механизмы оптического сглаживания неоднородностей. Добавление кластеров тяжелых элементов приводит к существенной доли лазерного излучения, преобразованного в рентгеновское, которое выносится вовне из мишени, при этом, естественно, уменьшается интенсивность газодинамического движения. Результаты показаний РЭОК можно увидеть на рис.18, и соответствующие результаты моделирования показаний РЭОК и баланс энергии – на рис.19.
База данных оптических коэффициентов DESOPLA содержит коэффициенты, равновесные по ионному составу, а также рассчитанные при учете неравновесности ионного состава плазмы. Как было показано Н.Н.Калиткиным в 2008 году, эффекты неидеальности плазмы очень слабо влияют на термодинамические функции, но оказывают влияние на ионный состав. Сравнение результатов при одной и той же энергии выстрелов с равновесными по ионному составу коэффициентами поглощения из базы данных DESOPLA и с неравновесными (рис.19), показывает, что при использовании неравновесных коэффициентов в расчете значительно увеличивается вынос энергии вовне и замедляется скорость рентгеновского фронта. Таким образом, использование неравновесных по ионному составу коэффициентов в большей мере отвечает физической ситуации в эксперименте.
Сводные данные по вложенным в пористые мишени энергиям и скоростям рентгеновского и гидротеплового фронтов, наблюдаемых в экспериментах и полученных в расчетах, представлены в Таблице 5.
а) выстрел № 28211, Elas= 158,6 Дж б) выстрел № 28220, Elas= 155 Дж Рис. 18. Показания РЭОК для пены TAC+9.9%Cu (9.1 мг/см3, 3 гармоника лазера, фольга 5мкм) в эксперименте.
a) Eabs = 130Дж;
1нс. Равновесные b) Eabs = 130Дж коэффициенты поглощения c) Eabs = 130 Дж;
1нс. d) Eabs = 130Дж Неравновесные к-ты поглощения e) Eabs = 18 Дж;
1нс. f) Eabs = 18 Дж Равновесные к-ты поглощения.
Рис. 19. Моделирование показаний РЭОК (слева) и энергетический баланс системы (справа) для пены TAC+9.9%Cu (9.1 мг/см3, 3 гармоника лазера, фольга 5мкм) в расчетах при разных энергиях выстрелов.
Таблица 5.
Сравнение экспериментальных и численных результатов по скорости рентгеновского фронта в зависимости от выбранной поглощенной энергии.
Vx (107 см/с) Vht (107 см/с) № выстрела E, Дж Свойства мишени ср. данные 7.7 2. 28205 170 6.6 2. 28207 157.5 11. 1. 28218 94 10. 2. 9,1 мг/см3, 28232 163 9. 3. моделирование 130 13.1 5. моделирование 50 7.2 3. моделирование 6 3–5.4 2. 28268 173 2. 28256 174 3. 4,5 мг/см3, моделирование 50 4. моделирование 10 ср. данные 7.5 2. 28211 159 8 2. 28220 155 6.4 TAC+ Cu 9,1 мг/см3, 3 моделирование 130 Eq 12.4 5. моделирование 130 NE 11. 4. моделирование 18 3.6 1. Таблица 5 показывает, что скорости рентгеновского и гидротеплового фронтов, полученные в экспериментах, измеряются неустойчиво. Сравнение теоретических и экспериментальных результатов позволяет сделать следующие выводы. 1) Перенос излучения определяет перераспределение энергии в системе. 2) Использование неравновесных по ионному составу коэффициентов поглощения в большей мере отвечает результатам экспериментов. 3) Величина поглощенной лазерной энергии составляет 30-40% от заявленной в эксперименте.
Для примера приведем сравнение для аналогичных экспериментов на установке LIL.
Поток лазерной энергии на поверхности мишени того же порядка, что и на установке PALS.
Существенное отличие заключается в длительности импульса (около 3 нс) и размере пятна на поверхности мишени (радиус 500 мкм). Приведем данные по полной лазерной энергии и пропущенной энергии, полученные в эксперименте и в расчете (рис. 20 и рис. 21).
Рис.20. Временные зависимости падающей, Рис.21. Временные зависимости падающей и пропущенной и отраженной лазерной пропущенной лазерной мощности при мощности для эксперимента ILP3. моделировании эксперимента ILP3.
Плотность TMPTA однородна.
Оценки величины компонентов баланса энергии, сделанные авторами экспериментов на LIL, практически точно совпадают с результатами моделирования (см. рис. 22), кроме одного пункта: авторы оценивали потери на рентгеновское излучение нагретой плазмы в Дж и признали, что в своих оценках потеряли по крайней мере порядок величины. При математическом моделировании баланс энергий сходится и представлен на рис. 22.
Нетрудно видеть, что потери на излучение составляют значительно большую величину, чем в физических оценках, – 1200-1500 Дж.
Рис. 22. Баланс энергий для моделирования эксперимента ILP3 на установке LIL.
Результаты Главы 4 опубликованы в работах [18-25].
В Главе 5 описан метод пересчета усредненных по спектру сечений для динамического моделирования саморегулируемых нейтронно-ядерных режимов. Методы эффективного понижения размерности уравнения переноса на основе квазидиффузионного подхода позволяют получить эффективную одногрупповую систему уравнений для полного скалярного и векторного потока нейтронов, которая может быть объединена в динамическом расчете с уравнениями выгорания, реакторной кинетики и управления.
Исследование саморегулируемых нейтронно-ядерных режимов (СНЯР) в активных зонах (АЗ) быстрых реакторов требуют динамического моделирования процессов, происходящих в АЗ. Эти процессы описываются уравнениями переноса нейтронов, выгорания и реакторной кинетики, а также введенного управления. Идея СНЯР первого рода (СНЯР-1) была предложена Л.П.Феоктистовым в 1988 году для гомогенной зоны. Она заключалась в том, что если в АЗ критическая концентрация плутония-239 ниже равновесной, то может возникнуть саморегулируемый режим, в котором концентрация плутония-239 со временем слегка растет, подтягиваясь к равновесной. В.Я.Гольдиным с сотрудниками в 1992г. были начаты работы по математическому моделированию СНЯР.
Была подтверждена возможность реализации СНЯР-1, и в 1995г. был предложен СНЯР-2, основанный на гетерогенности АЗ. Введение зон малого обогащения (в которых концентрация плутония ниже равновесной) и зон большого обогащения (в которых концентрация выше равновесной) в АЗ позволяет значительно улучшить параметры саморегулируемого режима по длительности кампании, равномерности энерговыделения во времени и пространстве. Предложенное управление, при котором выведение поглотителя (карбида бора) из АЗ используется только в самом начале кампании для вывода реактора на режим, а затем используется только тонкое управление введением соединений обедненного урана, позволяет осуществить работу АЗ реактора в СНЯР-2 в течение трех лет и более, что существенно превышает реакторную кампанию существующих и проектируемых быстрых реакторов, составляющую 140 суток. При этом (кроме первых 10-20 суток до установления режима) реактор работает без запаса реактивности.
Так как было совершенно очевидно, что создаваемая модель должна быть динамической и включать уравнения выгорания и кинетики, то первоначальная модель была одномерно-цилиндрической, позволявшей хорошо представить главную физическую суть явлений. Надо четко представлять себе, что и в настоящее время полный нестационарный многогрупповой расчет уравнения переноса нейтронов совместно с уравнениями выгорания и кинетики с учетом запаздывающих нейтронов (и других процессов) в многомерных геометриях является чрезвычайно трудоемким делом, возможным в единичных ситуациях на высокопроизводительных машинах. Трехмерные расчеты ядерных реакторов чаще всего сводятся к решению стационарных уравнений диффузии, реже – переноса. Метод квазидиффузии позволяет эффективно понижать размерность задачи, сводя ее к нестационарной усредненной одногрупповой, численное решение для которой совместно с уравнениями выгорания и кинетики ищется на каждом шаге по времени. При этом многогрупповая система уравнений переноса нейтронов с квазидиффузией в квазистационарном приближении время от времени пересчитывается, так что усредненные микросечения реакций со временем также изменяются в соответствии с изменением спектра.
Улучшение математической модели велось постепенно: сначала были введены в расчет двумерные утечки через торцы реактора на основе расчета двумерного одногруппового уравнения диффузии, а потом и одногруппового уравнения переноса с квазидиффузией, усредненные коэффициенты для которого были получены в ходе одномерного расчета. В диссертации описывается переход к двумерной многогрупповой задаче решения уравнения переноса и получение эффективной одногрупповой системы уравнений квазидиффузии для дальнейшего использования в нестационарном двумерном расчете совместно с уравнениями выгорания и кинетики. Заметим, что в предлагаемом СНЯР-2 для реактора реального типа отпадает необходимость в тяжелой системе компенсирующих стержней, управление реактором значительно облегчается, поэтому геометрия реактора лишается сильной асимметрии и практически близка к двумерной цилиндрической геометрии.
Расчеты в настоящей работе выполнены для стандартного 26-группового приближения библиотеки БНАБ-93.
Поскольку саморегулируемый нейтронно-ядерный режим поддерживает в реакторе p с очень малой квазистационарное распределение скалярного потока нейтронов постоянной времени, будем решать стационарное уравнение переноса нейтронов, которое в r-z-геометрии имеет вид:
p p p p N w + ( N tp + p ) p = ( )k d + z + r + k p r s k z r v k =1 (15) P 1 N (1 ) + 4 c + +Q k k k p k p p f f k gg dg k =1 g = Здесь угловые переменные азимутальный угол [ 0, ], полярный угол [ 0, 2 ], d = sin d d, z = cos, r = sin cos, = sin sin, (16) Первый член справа в (15) отвечает рассеянию всех видов из высокоэнергетичных групп и внутри рассматриваемой группы;
для быстрых реакторов рассеяния с увеличением энергии нет. Остальные два члена описывают размножение нейтронов, включая учет запаздывающих нейтронов;
Qp – возможный внешний источник нейтронов.
L N = N l l, L – количество Для всех процессов предполагается где l = рассматриваемых элементов в цепочке превращений, включая все изотопы, Nl – концентрация изотопа номера l.
В соответствии с общепринятыми обозначениями:
• t= s + f + n + n2n + n3n – полное сечение всех процессов столкновения нейтронов в группе p (p=1,…,P), включая:
• s= s(e) + s(in) – сечение рассеяния (упругое и неупругое);
• n2n –сечение поглощения нейтрона с рождением двух нейтронов;
• s – суммарное сечение рассеяния и реакций n2n, n3n;
• wkp ( ) – индикатриса рассеяния;
• f – сечение деления;
• n – сечение поглощения нейтрона с излучением кванта;
• k – доля запаздывающих нейтронов в группе k;
• g/ – доля группы g в запаздывающих нейтронах;
• kp – доля нейтронов деления группы k, попадающих в группу p;
• dgp – доля запаздывающих нейтронов группы g, попадающих в группу p;
• cg – концентрация предшественников запаздывающих нейтронов;
• g = ln 2 / Tg, Tg – период полураспада предшественников запаздывающих нейтронов группы g.
g =.
Условие нормировки:
g = В диссертации описаны квазидиффузионные преобразования членов рассеяния и деления в правой части уравнения переноса, позволяющие привести это уравнение к виду p p 1p p N s {w(0) k + 3w(1) k ( zWz + rWr ) + + ( N tp + p ) p = z + r + k p k p k k r 4 k = z r v } + 5 / 2 w(2) k ( 3( 2 Drr + 2 r z Drz + 2 Dzz + D ) 1) k + p k k k 2 k r z + N ( f f f f ) p + N f f ( dp f ( p ) ( )) + Q p g p f ( p ) ( ) = dg. (17) g =1 + g В этой форме уравнений важно, что члены рассеяния и деления приводятся к виду, содержащему полные и групповые скалярный и векторный потоки, а также ряд коэффициентов, усредняемых по спектру решения (члены с чертой сверху). При этом в соответствии с методом квазидиффузии групповые скалярные и векторный потоки нейтронов находятся из независимой системы многогрупповых уравнений квазидиффузии.
Эта система для задачи на собственные значения имеет вид (18)-(20), а для пересчета сечений при известной постоянной времени имеет вид (21)-(23):
Wz p 1 (rWr p ) + N ( tp sp w(0) p ) p = + p z r r (18) {( )} p = N w + (1 ) N f f + f f ( ), k p k p p p s (0) k d k = ( Dzz p ) 1 (rDzr p ) p p p + N ( t s w(1) p )Wz = N sk w(1) kWzk, + p pp p p (19) z r r k = ( Dzr p ) ( Drr p ) 2 Drr + Dzz 1 p p p p p p + N ( tp sp w(1) p )Wr p = N sk w(1) kWrk.
+ + p p (20) z r r k = Wz p 1 (rWr p ) + ( N ( tp sp w(0) p ) + p ) p = + p z r r v (21) p = N w + N ( f f f f ) + N f f ( f ( )), k p k p p ( p) s (0) k d k = ( Dzz p ) 1 (rDzr p ) p p p p Wz = N s w(1) kWz, + N ( tp sp w(1) p ) + p + p kp k (22) z r r v k = ( Dzr p ) ( Drr p ) 2 Drr + Dzz 1 p p p p p p p Wr = N s w(1) kWr.
+ N ( tp sp w(1) p ) + p + + p kp k (23) z r r v k = В методе квазидиффузии (при малой анизотропии рассеяния) главные члены рассеяния внутри самой p-ой группы переносятся налево для нулевого и первого моментов разложения индикатрисы рассеяния по полиномам Лежандра, т.е. учитывается их главная часть.
Решение для скалярного и векторного потока зависят от второго (и, возможно, более высоких) моментов индикатрисы рассеяния только через коэффициенты квазидиффузии, поэтому итерации по рассеянию объединяются с итерациями члена деления. В правой части уравнений (18)-(20) или (21)-(23) первая сумма отвечает рассеянию из энергетически более высоких групп, и к моменту расчета данной группы p эти члены на данной итерации источника уже известны. Этим расчет реактора на быстрых нейтронах выгодно отличается от аналогичного расчета реакторов на тепловых нейтронах, где существует рассеяние из низкоэнергетических групп в группы с большей энергией. Остальные два члена деления в (18) или (21) не могут быть определены до решения системы (18)-(20) или (21)-(23) для всех групп, поэтому они берутся из решения усредненных по энергии уравнений квазидиффузии с предыдущей глобальной итерации.
Усреднение в одну группу (21)-(23) производится по аналогии с методом, который был разработан для задач переноса излучения. Полученная система эффективных одногрупповых уравнений квазидиффузии может быть записана не в квазистационарном приближении, а в нестационарном виде, пригодном для динамического расчета (усреднение уравнений для поиска собственного значения аналогичны) 1 Wz 1 (rWr ) + N ( f + n n 2 n 2 n 3n ) = N ( f f f f f ( )) + + (24) v t z r r ( Dzz ) 1 (rDzr ) + N ( t s w(1) )Wz = d z, + (25) z r r ( Dzr ) ( Drr ) 2 Drr + Dzz + N ( t s w(1) )Wr = d r.
+ + (26) z r r Правые части dz и dr включают в себя как просуммированные по всем группам правые части уравнений (19)-(20), так и добавки, связанные с возможным знакопеременным усреднением в левых частях уравнений. Для ряда задач эти добавки не очень существенно влияют на расчет критических параметров, но с увеличением геометрических размеров АЗ быстрых реакторов влияние этих поправок становится существеннее.
Таким образом, при нахождении собственного значения вместо двойного цикла итераций (по рассеянию и делению и для нахождения или Keff) метод квазидиффузии позволяет использовать один, причем быстро сходящийся. Проведенные расчеты активных зон быстрых реакторов типа БН-800 и БОР-60, способных работать в СНЯР-2, потребовали 8-20 итераций для нахождения Keff. В обычно используемых методах источника суммарное число итераций по рассеянию и делению (т.е. количество решений многогруппового уравнения переноса для определения Keff) значительно больше. Нахождение критических параметров начальной сборки делает необходимым еще один итерационный процесс приведения Keff к единичному значению методом секущих.
На основании предложенных методик были рассчитаны критические параметры реакторов БОР-60 и реактора типа БН-800, которые способны работать в СНЯР-2.
Оптимизация режима СНЯР-2 по длительности кампании, равномерности энерговыделения во времени и пространстве требует большого количества динамических расчетов. Для их проведения используется более дешевая полуторамерная модель, использующая двумерные утечки через торцы реактора. В диссертации предложено ее улучшение введением не только средних утечек, но и утечек в каждой группе.
Результаты Главы 5 опубликованы в работах [26-31].
Основные положения, выносимые на защиту.
1. Созданы эффективные методики и комплексы программ решения многогруппового уравнения переноса совместно с квазидиффузией для решения задач переноса излучения в сплошной среде. Предложенные методы обладают повышенными свойствами монотонности и учитывают особенности решения. Методы эффективного понижения размерности уравнения переноса позволили создать экономичную и точную методику, учитывающую взаимное влияние переноса фотонов и газодинамических процессов в системе.
2. Предложен метод учета сильной анизотропии рассеяния в уравнении переноса, значительно сокращающий количество итераций по рассеянию. Решен ряд задач атмосферной радиации с рассеянием на аэрозолях и облаках, обладающих особенностью преимущественного рассеяния вперед. Применение предложенного метода учета анизотропии рассеяния совместно с методом лебеговского усреднения по частоте (А.В.Шильков) позволило получить прецизионные результаты, не имевшие аналогов в мире, для задачи об энергетическом балансе атмосферы Земли.
3. На основании разработанных автором методик расчета переноса излучения и известных газодинамических методик создан программный комплекс LATRANT для моделирования задач радиационной газовой динамики в r-z-геометрии при существенной роли собственного излучения плазмы. Полномасштабное моделирование задач УТС позволило объяснить экспериментальные результаты, полученные на установках PALS и LIL.
4. Создан эффективный метод и комплекс программ расчета многогрупповой системы уравнений переноса нейтронов с квазидиффузией в двумерной r-z геометрии, значительно сокращающий число итераций по рассеянию и делению, применяемый для проведения поисковых работ по оптимизации активных зон быстрых реакторов нового типа, предложенных и разрабатываемых в ИММ РАН, которые обладают повышенными свойствами безопасности и экономичности.
ПУБЛИКАЦИИ АВТОРА ПО ТЕМЕ ДИССЕРТАЦИИ 1. Е.Н. Аристова, Д.Ф. Байдин, В.Я. Гольдин. Два варианта экономичного метода решения уравнения переноса в r-z геометрии на основе перехода к переменным Владимирова // Математическое моделирование, 2006, т. 18, № 7, стр.43-52.
2. D.F. Baydin, E.N. Aristova, V.Ya. Gol’din. Comparison of the efficiency of the transport equation calculation methods in characteristics variables // Transport Theory and Statistical Physics, 2008, v.37, № 2&4, pp. 286-306.
3. Е.Н. Аристова, А.В. Колпаков. Комбинированная разностная схема для аппроксимации эллиптического оператора на косоугольной ячейке // Математическое моделирование, 1991, т.3, №4, стр.93-102.
4. Д.Ю. Анистратов, Е.Н. Аристова, В.Я. Гольдин. Нелинейный метод решения задач переноса излучения в среде // Математическое моделирование, 1996, т.8, №12, стр.3 29.
5. Е.Н. Аристова, В.Я. Гольдин, А.В. Колпаков. Методика расчета переноса излучения в теле вращения // Математическое моделирование, 1997, т.9, №3,с.91-108.
6. Е.Н. Аристова, В.Я. Гольдин, А.В. Колпаков. Перенос излучения через кольцевую щель в теле вращения // Математическое моделирование, 1997, т.9, №4, с.3-10.
7. E.N. Aristova, V.Ya. Gol’din, A.V. Kolpakov. Multidimensional Calculations of Radiation Transport by Nonlinear Quasi-Diffusional Method. Proceedings of the Joint Intern.Conference M&C;’99 on Mathematics and Computation, Reactor Physics and Environmental Analysis in Nuclear Applications, Published by Senda Editorial, S.A. Isla de Saipan, 47, 28035 Madrid, Spain, v.1, p.667-676.
8. Е.Н. Аристова, В.Я. Гольдин. Эффективное понижение размерности уравнения переноса. Энциклопедия низкотемпературной плазмы, 2000, Вводный том, т. 1, с. 462-471.