Математическое моделирование переноса излучения и переноса нейтронов с учетом процессов в сплошных средах

Размещено на http://www.allbest.ru/

На правах рукописи

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

Автореферат

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

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

Аристова Е.Н.

Москва - 2009 год

Работа выполнена в Институте математического моделирования РАН

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

член-корреспондент РАН, профессор Суржиков Сергей Тимофеевич

доктор физико-математических наук, профессор Трощиев Виталий Ефимович

доктор физико-математических наук, профессор Шагалиев Рашит Мирзагалиевич

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

Физический институт им. П.Н.Лебедева РАН

Защита состоится « 9 » апреля 2009г. в 1100 часов на заседании диссертационного совета Д.002.058.01 при Институте математического моделирования РАН по адресу: Москва, 125047, Миусская пл, д. 4а.

С диссертацией можно ознакомиться в библиотеке ИММ РАН.

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

Ученый секретарь диссертационного совета Д.002.058.01, доктор физико-математических наук Змитренко Н.В.

1. Общая характеристика работы

Актуальность темы

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

Трудности решения уравнения переноса помимо большой размерности связаны с несколькими факторами:

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

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

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

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

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

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

Цель и задачи исследования

Настоящая работа посвящена разработке эффективных численных методов решения уравнения переноса в рамках квазидиффузионного подхода. Метод квазидиффузии (В.Я. Гольдин. Квазидиффузионный метод решения кинетического уравнения // ЖВМ и МФ, 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 страницах.

2. Основное содержание работы

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

Метод квазидиффузии заключается в постепенном понижении размерности задачи так, что на последнем этапе становится возможным эффективное объединение следствий уравнения переноса с уравнениями сплошной среды, описывающими другие физические процессы, например, с уравнениями газовой динамики для излучения или реакторной кинетики и выгорания для переноса нейтронов. Усреднением многогрупповых уравнений переноса по углам получается многогрупповая система уравнений квазидиффузии, замкнутая при помощи дробно-линейных функционалов, вычисляемых по решению уравнения переноса. Усреднение полученной многогрупповой системы уравнений квазидиффузии по энергии приводит к эффективной одногрупповой системе уравнений квазидиффузии относительно полных плотности и потока излучения, которая может быть объединена с другими уравнениями сплошной среды, что решает главную задачу - учета нелокального и нелинейного взаимодействия излучения и вещества. Внешнее усложнение процедуры решения уравнения переноса позволяет разрешить ряд трудностей. Во-первых, использование многогрупповой системы уравнений квазидиффузии автоматически означает консервативность схемы. Во-вторых, выделение главной части в члене рассеяния внутри группы позволяет перенести ее налево для уравнений квазидиффузии, так что влияние остаточных членов становится опосредованным через тензор квазидиффузии. Итерационный процесс решения уравнений с такой иерархической структурой быстро сходится при условии малости производных Фреше от введенных дробно-линейных функционалов по решению. Таким образом, решение уравнения переноса сводится к ряду вычислительных задач: 1) решение уравнения переноса при известной правой части, усреднение полученного решения по углам и вычисление тензора квазидиффузии, 2) решение многогрупповой системы уравнений квазидиффузии и усреднение этой системы по энергии в одногрупповую систему, 3) совместное решение полученной одногрупповой системы уравнений квазидиффузии с другими уравнениями сплошной среды (уравнениями энергии для задач радиационной газовой динамики или выгорания и реакторной кинетики). Диссертация посвящена методам решения поставленных вычислительных задач и моделированию физических задач на основе разработанных алгоритмов.

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

,

, , , , .

Переход к собственным характеристическим переменным позволяет записать уравнение, содержащее только две пространственные производные:

,

Рис.1С формальной математической точки зрения связь r и ц с переменными s и h осуществляется по формулам (рис.1)

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

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

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

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

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

Результаты Главы 1 опубликованы в работах [1-2].

Глава 2 посвящена созданию и тестированию всех методов и программ, необходимых при численном решении системы квазидиффузионных уравнений и объединении их с уравнением энергии для вещества, что включает в себя: 1) построение разностных схем для многогрупповой системы уравнений квазидиффузии, 2) разработку методов решения полученных разностных уравнений эллиптического типа, 3) усреднение уравнений в эффективную одногрупповую систему, 4) решение объединенной системы усредненных уравнений квазидиффузии с уравнением энергии для вещества.

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

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

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

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

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

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

.

Данная задача характерна тем, что фронт излучения в основной части трубы перпендикулярен контактной границе. На контактной границе происходит изменение направления движения волны излучения от распространения вдоль оси z для внутренности трубы к распространению по радиусу для стенок трубы. Это приводит к большим значениям недиагонального компонента тензора квазидиффузии около контактной границы (рис. 3 б), которыми нельзя пренебречь для построения аналога монотонной схемы. Использование аналога монотонной схемы с поворотом координат в плоскости r-z позволяет провести расчет этой задачи.

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

Температурные поля для пяти моментов времени приведены на рис.5.

Результаты Главы 2 опубликованы в работах [3-11].

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

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

Сходимость метода тестировалась на задаче рассеяния в однородном плоском слое 0?z?10, граничащем с вакуумом с двух сторон, с изотропным источником q=1, полным сечением у=1 и сечением рассеяния уs=с. Для индикатрисы рассеяния Хеньи-Гринстейна

,

была исследована зависимость числа итераций по рассеянию от параметров жесткости задачи g и c. Здесь g - параметр, отвечающий за особенность преимущественного рассеяния вперед. Чем ближе величина g к единице, тем ближе особенность к д-особенности при м>1. В обычных методах итераций источника наихудшая сходимость имеет место при отсутствии поглощения. Количество итераций по рассеянию в предложенном методе (последние четыре колонки) в сравнении с рядом других методов из работы А.В.Волощенко представлено в Табл. 1.

Наибольшее число итераций по рассеянию возникает для задачи средней жесткости (параметр g=0.99) при отсутствии выделения д-особенности. В остальных случаях число итераций невелико и колеблется от 9 до 15. Ужесточение критерия сходимости ведет к двукратному росту числа итераций при замене е на е2.

Таблица 1. Число итераций при критерии сходимости при изменении параметров g и c.

g	c	Число итераций
		P1SAa	N/Ab	P1SA-дc
					0.95	0.99	0.999	1.
0.9	0.9 0.99	20-15 66-23	66-45 392-89		9 12	11 17		12 19
0.99	0.9 0.99	28 75	48 120	13-21			12 15	21 21
0.999	0.99			5-7			12

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

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

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

Рис.6а. Индикатриса рассеяния для модели водяного облака с широким распределением частиц по размерам, обладающая сильной особеннос-тью преимущественного рассеяния вперед.. Рис.6б. Выделение -особенности в индикатрисе рассеяния. Пунктиром нарисована регулярная компонента wreg().

Стационарное уравнение переноса в приближении плоской атмосферы для света заданной частоты имеет вид:

,

где - коэффициент поглощения, а - коэффициент рассеяния.

В предлагаемом методе интенсивность излучения I представляется в виде суммы трех компонент:

1) I0 - нерассеянное солнечное излучение, в том числе зеркально отраженное от поверхности земли,

2) I1 - излучение, рассеянное на малые углы, в том числе с многократным рассеянием,

3) I2 - все оставшееся излучение.

На верхней границе атмосферы задается падающее солнечное излучение в виде -функции угла: I=q·(-0). На нижней границе атмосферы, т.е. на поверхности земли с альбедо r, задается условие отражения, при этом предполагается, что доля излучения отражается зеркально, а доля 1- - диффузно. Метод позволяет учитывать и более сложные варианты граничных условий отражения, но за неимением соответствующих данных пришлось ограничиться комбинацией зеркального и диффузного отражения.

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

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

Уравнение переноса для компоненты I2 содержит решения для I0 и I1 . Для его решения используется метод квазидиффузии, модифицированный для учета сильной анизотропии рассеяния.

Отметим характерные черты предложенного метода.

Вычисление интеграла рассеяния.

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

Для корректного учета рассеяния на малые углы индикатриса рассеяния делится на гладкую регулярную часть с областью определения [-1,1], и сингулярную, с малым носителем в окрестности =1: w()=wreg()+wsing() (рис.6б). Главная часть сингулярной компоненты индикатрисы рассеяния может быть представлена в виде -функции, подстановка которой в уравнение переноса приводит к эффективному уменьшению коэффициента рассеяния. При этом в правой части остается поправочный член рассеяния на малые углы, который можно не учитывать в тех областях фазового пространства, где рассеяние на малые углы несущественно (вдали от сосредоточенного источника).

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

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

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

Если выбран метод вычисления интеграла рассеяния от остатков индикатрисы с использованием замены переменных, сетка интегрирования по продуцируется сеткой, на которой задана индикатриса рассеяния (см. пп. а)), в силу однозначности связи и при фиксированных ,/ ,[0,].

Член рассеяния на малые углы после выделения главного члена wsing в предположении дифференцируемости интенсивности излучения может быть приведен к Фоккер-Планковскому виду, т.е. к диффузии в фазовом пространстве. Однако мы не используем это представление, поскольку с нашей точки зрения предположение о дифференцируемости функции распределения является очень сильным предположением, которое вряд ли справедливо в окрестности =0 и =0.

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

Численная схема как для уравнения переноса, так и для уравнений квазидиффузии выписывается в предположении постоянства коэффициентов поглощения и рассеяния на каждом расчетном интервале. Нерассеянное солнечное излучение входит в виде экспоненциального источника в уравнения для компонент I1, I2, кроме того, в уравнения для I2 входит решение для компоненты I1. Численная схема учитывает экпоненциальную зависимость I0, а зависимость всех остальных членов по пространству предполагается кусочно-линейной.

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

Была решена задача о рассеянии солнечного излучения для реальной модели атмосферы, состоящей из воздуха с рэлеевским законом рассеяния и ряда модельных аэрозолей (фонового стратосферного, континентального, городского, морского - рис. 10), при отсутствии поглощения. Рассеяние на аэрозолях также обладает особенностью преимущественного рассеяния вперед, соответствующие индикатрисы рассеяния были рассчитаны в книге (М.Я.Маров, В.П.Шари, Л.Д.Ломакина. Оптические характеристики модельных аэрозолей атмосферы земли. ИПМ им. М.В.Келдыша АН СССР, М.: 1989, 230с) и представлены на рис.7. Из предложенных для сравнения вариантов мы выбрали расчет для длины волны 0.55 при нормальном падении солнечного излучения. Распределение аэрозолей соответствует профилю II рабочей группы «Стандартная радиационная атмосфера» (Перенос радиации в рассеивающих и поглощающих атмосферах. Стандартные методы расчета. Под ред. Ж.Ленобль. - Л.: Гидрометеоиздат, 1990, 230с.) и представлено в Табл. 2. Все остальные данные также соответствуют этому варианту. Расчет велся в физических переменных в слое 104км0км.

Таблица 2. Распределение аэрозолей и оптические толщины для реальной модели атмосферы.

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



Рис.7. Задача IV. Индикатрисы рассеяния континентального, морского, городского и стратосферного аэрозолей с увеличенным масштабом углов преимущественного рассеяния вперед

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

Косинус угла падения солнечного излучения	Результат предложенного метода Поток постоянен по высоте	Результаты рабочей группы
		z=30км	z=2км	z=0км
м0=1.	2.575	2.863	2.852	1.785
м0=0.75	1.759	2.070	2.059	1.097
м0=0.5	0.997	1.292	1.282	0.517

А.В.Шильковым с соавторами была разработана система кодов атмосферной радиации ATRAD, включающая восстановление коэффициентов поглощения атмосферными газами по параметрам линий из базы данных HITRAN и подготовку коэффициентов поглощения и рассеяния, усредненных по Лебегу. Вид уравнения переноса сохраняется при использовании лебеговских групп. Точность метода лебеговского усреднения на 400 эффективных группах совпадает с точностью методов типа `line-by-line' при сокращении на четыре порядка количества используемых энергетических точек. Разработанная методика учета анизотропии рассеяния была включена в систему кодов атмосферной радиации ATRAD.

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

Таблица 4. Полное поглощение всеми атмосферными газами в линиях с учетом и без учета молекулярного рассеяния (Вт/м2).

	Альбедо поверхности 0.2	Альбедо поверхности 0.8
	И=30є	И=75є	И=30є	И=75є
Поглощение	178.2	73.8	200.4	81.7
Поглощение и рассеяние	175.6	71.0	197.5	79.0

Изменение потока по высоте для случая поглощения всеми атмосферными газами без учета и с учетом рассеяния приведено на Рис. 8. Результатов `line-by-line' расчетов для рассеивающих атмосфер крайне мало.

На рис. 9 приведена зависимость скорости радиационного выхолаживания от давления, полученная методами `line-by-line' разными авторами для чисто поглощающей атмосферы и по расчетам в системе ATRAD.

Таким образом, можно заключить, что метод лебеговского усреднения дает точность расчета методов `line-by-line' при сокращении объема вычислений в 10000 раз. Предложенный метод учета анизотропии рассеяния хорошо работает как в лебеговских обобщенных, так и в физических переменных. При объединении метода лебеговского усреднения с методом учета сильной анизотропии рассеяния удается получить прецизионные результаты расчетов по тепловому балансу атмосферы Земли с учетом рассеяния, которые практически недоступны для других методов расчета.

Результаты Главы 3 опубликованы в работах [12-17].

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

Глава 4 посвящена описанию программного комплекса LATRANT, созданного для математического моделирования задач инерциального термоядерного синтеза. Программный комплекс учитывает перенос излучения в многогрупповом приближении и движение газа по улучшенной лагранжевой методике в двухтемпературном приближении в переменных r-z. Программный комплекс был создан на базе газодинамического комплекса ATLANT, созданного А.Б.Искаковым, и программного комплекса LATRA совместного решения уравнений переноса излучения и энергии вещества, созданного автором диссертации.



Рис.9. Скорость радиационного выхолаживания для стандартной летней атмосферы средних широт.

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

,

,

.

Здесь:



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

отвечает за кинетику ионизации, в этой формуле Z - заряд иона, I(Z) - энергия ионизации, - скорость ионизации;

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

планковская равновесная плотность излучения

(h - постоянная Планка, k - постоянная Больцмана,

у - постоянная Стефана-Больцмана),

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

Система (1) - (7) замыкается уравнениями состояния:

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

Для спектрального описания задачи используется многогрупповое приближение, уравнение переноса в квазистационарном приближении в r-z геометрии имеет вид:

,

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

,

последний из коэффициентов используется в многогрупповых уравнениях квазидиффузии. Здесь - эффективная температура усреднения. Используемая в расчетах база данных DESOPLA (ФИ РАН) включает в себя таблицы групповых коэффициентов поглощения , , для значений плотностей 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 кэВ для наиболее часто используемых в экспериментах по УТС веществ. Под первоначально понималось максимальное значение температуры в расчете, которое должно было либо априорно оцениваться, либо определяться из предварительного расчета. В дальнейшем была сделана модификация использования всей базы данных по , используя оценку величины по положению максимума пришедшего излучения.

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

Расщепление по физическим процессам приводит к следующей схеме расчета уравнений на одном временном шаге:

Запоминание всех необходимых нестационарных величин с предыдущего временного шага;

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

Расчет вклада электронной и ионной теплопроводности в уравнения энергии;

Решение групповых уравнений переноса, их усреднение по направлениям полета фотонов в квазидиффузионные многогрупповые уравнения;

Суммирование и усреднение групповых уравнений квазидиффузии по энергии в эффективную одногрупповую систему уравнений квазидиффузии;

Совместное решение двух уравнений энергии для электронного и ионного компонента плазмы и эффективной одногрупповой системы уравнений квазидиффузии;

Расчет энергетического баланса.

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

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

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

...