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

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

Нифанова Александра Васильевна

ХАРАКТЕРИСТИЧЕСКИЕ Бп-МЕТОДЫ ДЛЯ КИНЕТИЧЕСКОГО УРАВНЕНИЯ ПЕРЕНОСА НЕЙТРОНОВ В СФЕРИЧЕСКИХ СИСТЕМАХ

Специальность 05 13 18 - математическое моделирование, численные методы и комплексы программ

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

Троицк-2008 003444а 16

003444916

Работа выполнена в ФГУП «Государственный научный центр Троицкий институт инновационных и термоядерных исследований»

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

Официальные оппоненты доктор физико-математических наук, профессор Головизнин Василий Михайлович, ИБРАЭ РАН

кандидат физико-математических наук, старший научный сотрудник Аристова Елена Николаевна, Институт математического моделирования РАН

Ведущая организация Институт теоретической и математической физики РФЯЦВНИИЭФ

Защита состоится « 24 » сентября 2008 г в (Ь часов Ь0 МКН

На заседании диссертационного совета Д212 130 09 в Московском инженерно-физическом институте поадресу 115409,г Москва, Каширское шоссе, д 31,тел 324-84-98,323-92-56

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

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

Автореферат разослан « Я _» (ОМАЛ 2008 г Ученый секретарь диссертационного совета

доктор физико-математических наук, профессор Леонов А С

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

Диссертация посвящена построению новых разностных схем с улучшенными свойствами монотонности на базе классического БЗп-метода и метода характеристических трубок Новые схемы предназначены для решения задач переноса частиц в системах со сферической симметрией

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

Для широкого круга задач нейтронно-ядерной физики со сферической симметрией процессов переноса и кинетики нейтронов основными математическими методами их численного решения являются $п-методы и методы характеристик Первые варианты этих методов были независимо сформулированы в конце 1940-х и начале 1950-х годов в работах по атомным проектам США и СССР Это Бп-метод Карлсона [1] и КН-схема Гольдина [2], метод прямого интегрирования Рихтмайера [3] и метод характеристик Владимирова [4] Несколько позже были предложены дискретный Эп-метод (ОЗп-метод) [5] и метод характеристических трубок (ХТ-метод) [6], которые представляют собой развитие и обобщение в определенных направлениях первоначальных методов

Бп-методы - это конечно-разностные аппроксимации кинетического интегро-дифференци-ального уравнения переноса частиц, рассматриваемого как уравнение в частных производных первого порядка Методы характеристик - это разностные или разностно-аналитические аппроксимации семейства обыкновенных интегро-дифференциальных уравнений (ОИДУ), записанных на характеристиках-траекториях движения частиц в сфере Отсюда проистекает существенное различие в математических свойствах этих методов и в классах задач, для которых их применение эффективно

Дискретный Бп-метод является наиболее простым и экономичным с точки зрения программной реализации и объема вычислений Метод использует прямоугольные сетки (Бп-сетки) и всюду, кроме окрестности центра сферы, имеет второй порядок аппроксимации и точности на гладких решениях [7] ОБп-метод консервативен относительно законов сохранения нейтронов и легко обобщается на многомерные геометрии, но имеет большой недостаток, - метод немонотонен Это может приводить к осцилляциям в сеточном решении или даже к появлению отрицательных значений скалярного потока, что существенно снижает точность расчетов при численном решении стационарных и нестационарных задач переноса в сложных гетерогенных средах В разное время были предложены алгоритмы монотонизации [9—13], а также различные модификации ОБп-мстода [14, 15] для повышения его точности Они достаточно эффективны при определенных условиях, но решают проблему немонотонности лишь частично, а их обобщения на многомерные уравнения очень сложны

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

ч

ОБп- и ХТ-методы - принципиально разные по своей сути, но их аппроксимационной основой являются сеточные уравнения баланса, записанные соответственно на Бп- и Т-сетках Это ап-проксимационное свойство в определенной степени их сближает и может быть основой для построения новых численных методов

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

Цель работы

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

Основные результаты работы

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

2 Математические понятия инварианта переноса и среднего расстояния, ранее введенные в методе характеристических трубок обобщены на сетки произвольной формы Установлена их связь с фазовым объемом сеточных ячеек в сфере В Эп-ячейке, трактуемой как характеристическая трубка, построено ОДУ баланса относительно полного потока и функция независимого источника с непрерывным изменением аргумента - расстояния от освещенных граней до неосвещенных

3. Для ОДУ в Бп-ячейке предложена экономичная монотонная 2-го порядка точности разностная схема Эта схема вместе с различными алгоритмами распределения полного потока по неосвещенным граням представляют собой новый численный метод - ОЗп-метод характеристических трубок (В8п1-метод), в котором полностью устранена причина немонотонности, обусловленная аппроксимацией столкновительных членов в уравнении переноса и кинетики нейтронов

Ц Написаны программы, реализующие разработанные методы Эффективность новых ОБт-схем подтверждена численными расчетами задач с независимыми источниками и задач на собственные значения Полное устранение одной причины немонотонности приводит к качественно новым численным результатам (квазимонотонные схемы)

Достоверность результатов

Достоверность полученных результатов подтверждается сравнительными численными исследованиями новых схем и классического ОБп-метода, выполненными для различных классов

задач, а также сопоставлением с результатами исследований Бп-методов, проведенными ранее другими авторами [8, 9, 15—17] Расщепление причин немонотонности и устранение одной из них обосновано аналитическими преобразованиями

Научная новизна

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

Практическая значимость работы

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

Апробация и публикации

Основные результаты работы докладывались и обсуждались на Научных сессиях МИФИ (2001, 2003, 2004, 2007), на семинаре «Нейтроника-2005» в Обнинске, на семинаре В Я Гольдина Института математического моделирования РАН, на семинаре В М Головизнина ИБРАЭ РАН

По теме диссертации опубликовано 8 печатных работ, из них две статьи в реферируемых журналах - «Доклады академии наук», «Математическое моделирование», два препринта, тезисы докладов на Научных сессиях МИФИ

Структура и объем работы

Диссертационная работа состоит из четырех глав, введения, заключения, приложения и списка литературы Материал диссертации изложен на 97 страницах, включает 12 рисунков, 3 таблицы и список литературы из 58 наименований

Краткое содержание работы

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

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

г 8г дц

1-м

2 \ N

+ a(r)N = S(r)y

or г д/л

¿2-1

Краевые условия имеют вид N(R,M) = W(М)> И— О dN(r,ft = -l)

(1)

О') О")

dr

+а(r)N(г,ц = —l) = S(r), N(r = 0, //) = const

В уравнениях (1), (Г) г, ц - независимые переменные г - расстояние от центра сферы до точки М, где находится нейтрон, ц = cos в, где в- угол в точке М между радиус-вектором Г и направлением полета нейтрона £, в сфере (см рис 1) v - скорость движения нейтронов, t - время Таким образом, переменные г, fi изменяются в области D = {О < г < R—1 <//</}, где R -внешний радиус сферы Функции Р(г), <j(г), S(г) заданы

Искомая функция N(r,ju) есть плотность нейтронов в фазовом пространстве т,ц, иначе говоря, , dW = r2drd/j есть число частиц в элементарном фазовом объеме dW

Рис 1 Изменение угла в на траектории нейтрона в сфере, ¡1 = cos в

Рис 2 Область D

В задачах с источником Qfr.fi) ф О либо у(ц)фО В задачах на расчет критических параметров £~)(г,р) = 0 и

Уравнения переноса в недивергентной форме (Г) рассматриваются также как однопарамет-рическое семейство обыкновенных интегро-дифференциальных уравнений на характеристиках-траекториях полета частиц в пространстве

ад ад ад г

Применяется итерационный метод решения изложенных задач

Во второй главе проводится расщепление ОЭп-метода на схему метода характеристических трубок относительно полного потока частиц и схему распределения полного потока по неосвещенным граням Бп-ячейки

Классический ОБп-метод записывается в виде уравнения баланса частиц (2) в прямоугольной Бп-ячейке (рис 3,4) и дополнительных аппроксимационных соотношений (3)

Я ' •г

-1 0 +1

Рис 3 Бп-сетка в области О

1Р+1.Ч+11

Рис 4 Ячейка ¿„

0>/ ИР+, - 4 лОМд+1 + - - (/ - Аг +

+ аР+1 ЫР+1АуР+1 = Яр+1 А УР+1 Ам . (2)

1 1 о

АУР+1=\('"I♦1 -^ р = Р -1 > ? = ОД,-,? -1,

и г Э

(3)

/¿<0: искомые значения л^« и л^; ц>0 \ искомые значения и .

Метод характеристических трубок (ХТ-метод) записывается в виде уравнения баланса частиц (4) в характеристической ячейке-трубке (рис. 5, 6) относительно полных потоков фр и на ее торцах:

ФР+1 - Ф р + ~ сгр+1 {ЫР + А Рй = А р,ч, (4)

2 2 2

где ФР+, = гр+, /V Д , ФР = гР ИР Мр&мр, = Я г2 йгй^..

Лря

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

Рис. 5. Характеристическая ячейка-трубка.

-1 0 +1 I1

Рис. 6. Сетка ХТ-метода.

Далее В8п-метод преобразуется в каждой Бп-ячейке в эквивалентную двухэтапную схему расщепления. Для этого в ОЗп-методе по аналогии с методом характеристических трубок (4) [6]

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

,Ф/ + Фг

Ф2~Ф, + аР+1\ 2 2+а1 =

р+1

(5)

где

Ау" = Мя+1 - Мч. М^ = -2 Ьд+1+»д} А г = гР+, ~ гР, Гр+1 = - (гР+] + гр), А = у {гр+1 ~ гР)&М

ШР

— = А^ - среднее расстояние проходимое частицами в ячейке, а - величина 2-го порядка

малости относительно шагов Дг, А/л

Ф; = Ф/м = ^-^ГрМгИ, -

<$>2 = Ф2РЙ = {1~ гРД А гКд+1 - г\ А

Ф; = Фум = (/'- + /Р Ич+1

' для ц <0,

для ц>0,

Ф; = Ф 2РЧ = (7 - М2я+)гР+1 А + гР„

(7 ~ А г- г р = М*)Гр+1 А г- Г2Р+, М<0

(7 ~ А г+ гр+, Мч+1 А ^ в (/ - М2)гРЛ Дг+ гРМ^М. М>0-

Делается вывод о том, что уравнение (5) представляет собой конечно-разностную аппроксимацию (немонотонную и неположительную) некоторого обыкновенного дифференциального уравнения (ОДУ) относительно функции полного потока ф(£)

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

Во второй главе также рассмотрена возможность расщепления балансных разностных схем для многомерных кинетических уравнений переноса частиц [18]

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

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

Из уравнения баланса следуют выражения для полных потоков ф, и ф2 (6), фазового объема области 5<о и среднего расстояния <5£, проходимого в ней частицами, а также выражения для столкновительного члена и для правой части получаемого ОДУ баланса в зависимости от непрерывно меняющегося аргумента - расстояния от освещенного торца ячейки вдоль средней характеристики в ней (7)

Рис 7 Характеристическая трубка в области ц > О

-а

с,

-в,

-в,

ф2 = \NdiL (/ _ м)г2)=М21 .

(6)

<0 1

Зй) [\ы{ц.г\2(1г(1ц \\5{р,г)гЧгс1ц

\\r4rdn \Wdrdn (7)

В результате предельного перехода (при уменьшении области) в каждой Эп-ячейке выводится ОДУ баланса частиц

^ + <7(ф + /р) = 5/. 0</;< Д£,

и ставится задача Коши с начальным условием ф(о) = ф,, где р - величина 2-го порядка малости по ширине трубки

Алгоритм построения ОДУ баланса применяется к Бп-ячейке, которая рассматривается как ячейка-трубка

Рис 8 Бп-ячейка, характеристики и промежуточные контуры

Для Бп-ячейки записываются выражения всех необходимых величин как в случае ячейки-трубки общего вида, и формулируется задача Коши для ОДУ относительно функции полного потока через контур (8) или относительно - функции среднего значения плотности частиц на контуре (9)

+ аФ = 1Б, Ф(о) = фр 0 < ^ < Д^

(8)

После того как решена задача построения ОДУ баланса в Бп-ячейке дается обзор подходящих для его решения монотонных и условно-монотонных схем 2-го порядка точности [19—23] и выбирается из них для проведения численных расчетов наиболее простая и экономичная

(10)

Однако схема (10) является условно-монотонной Для монотонизации схемы (10) предложен и применяется простой алгоритм нелинейного типа Решение по монотонной положительной схеме ОДУ баланса есть первый этап предлагаемого в диссертационной работе ОЗп-метода характеристических трубок (ОБш-метод)

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

,4+1 - БЗп-распределение

Таким образом построены новые ОЗп-методы, состоящие на первом этапе из монотонной 2-го порядка точности схемы (10) для ОДУ баланса, а на втором - из 1-го или 2-го способа распределения полного потока (ОЗш-методы)

2

- ст Зп-распределение

\ I I

NN

Рис 9 Построение ст Эп-распределения

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

Задачи с источником. Сравнительный численный анализ новых двухэтапных схем с классическим DSn-методом для задач с гладкими решениями (однородные шары), в которых DSn-метод не проявляет немонотонности, показал, что DSnt-схемы по качеству ему не уступают.

Для задач с гладкими решениями, в которых DSn-метод на крупной сетке дает сильную немонотонность сеточного решения, применение новых схем на той же сетке существенно уменьшает немонотонность и обеспечивает практически 2-й порядок точности численного решения. Рассматривается задача с параметрами: R=20; a(r)=10; (3(r)=0; Q(r)=l. DSn-метод дает сильно немонотонное решение во всей области. DSnt-метод с DSn-распределением дает немонотонное решение, но немонотонность существенно уменьшается, это объясняется тем, что одна причина немонотонности устранена, и осредненный полный поток на неосвещенных гранях получен положительным, немонотонность остается только из-за способа распределения этого потока. DSnt-метод с aSn-распределением дает гладкое решение, полностью соответствующее физическому процессу. На

рис. 10 представлен скалярный поток п (/>) = X N р A/i , решение получено на равномерной сет-

я

ке (10 интервалов по переменной г, 8 интервалов по переменной ц). Обозначения на рисунке: DSN - Dsn-метод, CHDSN и СН_СН - новые DSnt-схемы, CHDSN - Dsn-распределение, СН_СН -aSn-распределение.

п» 0.15

0.10

о.оо

-DSN CH_DSN - сн сн

5.00 10.00 15.00 20.00

Рис. 10. Однородная сфера: R = 20; ст(г) = 10; J3(r) = 0; Q(r) = 1.

Для сложных многослойных задач новые схемы дают сеточное решение близкое к точному. В качестве примера приводится расчет тестовой задачи Рида.

Таблица 1. ТЕСТОВАЯ ЗАДАЧА РИДА

0 < г < 2 2 < г < 3 3 < г < 5 5 < г < 6 6 < г< 8

р 0,0 0,0 0,0 0,9 0,9

<т 50,0 50 0,0 1,0 1,0

0. 50,0 0,0 0,0 0,7 0,0

Рис. 11. Задача Рида, три метода решения.

Рис. 12. Задача Рида, ОБШ-метод с а5п-распределением.

Расчеты показывают, что устранение одной причины немонотонности приводит на практике к существенному улучшению численных результатов. Можно заключить, что двухэтапные ОБЩ-схемы являются квазимонотонными.

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

Рассматривается задача из работы [16] на определение наименьшего положительного числа а (а^), при котором однородная задача

ч дЫ 1-а дЫ 1Г 1 , . , ,

Ш\г,ц)= /2— +----+ аИ

дг г дц 2

N(l,ju) = 0, M<0 h = 1724

имеет ненулевое решение Расчет о^, был произведен DSnt-методом с ffSn-распределением на равномерной сетке (10 интервалов по г, 10 интервалов по ц), а*р = 1,280 Этот результат можно считать хорошим, в работе [16] приводится значение = 1,279, полученное методом характеристик и £Цф = 1,280, полученное методом Бубнова-Галеркина.

Также представлена задача определения критического параметра X слоистой системы, приведенная в работе[17]

( 1 Л

у дг дц г \ v)

V У

S(r) = Mllno(r)+£llli n°(r) = '\N(r,pi)dfd 2 2-i

6(r) = 0, N(5,0,¿i) = 0 для Ц^О

Сферически-симметричная система состоит из трех областей, заполненных средой с различными свойствами, параметры приведены в таблице 2

Таблица 2 ТРЕХОЕЛАСТНАЯСИСТЕМА v=100, Q(r) =0

0 < г < 1,0 1,0 < г < 3,0 3,0 < г < 5,0

а 1, 0,1 0,01

Р 2 0,1 0,01

Для численного расчета выбирались начальные сетки (такие же как в работе [17]) По переменной г 12 интервалов г0= 0,0, Г1=0,1, г2= 0,4, г3= 0,7, Г4=0,9, Г5=1,0, г6=1,2, г7 = 2,0, г8= 2,8, г9= 3,0, г10= 3,6, г„ = 4,8, г12=5,0

По переменной ц 10 интервалов До= —1,0, д( = -0,9, Цг=-0,7, Дз=-0,4, /¿4 =-0,1, ¡15 = 0,0, щ= 0,1, 0,4, д8= 0,7, ¡1, = 0,9, щ0= 1,0

Расчеты проводились также на пропорционально измельченных сетках в 2 и в 4 раза, значения X приведены в таблице 4 (Хэ - экстраполированное значение), результаты для непрерывного Зп-метода (кроме Хэ) взяты для сравнения из работы [17]

Таблица 3 РАСЧЕТ КРИТИЧЕСКОГО ПАРАМЕТРА \ПРИ ИЗМЕЛЬЧЕНИИ СЕТКИ

Число интервалов по г и по д 12 по г, lOno/i 24 по г, 20 по /1 48 по г, 40 по ¡1 Хэ

DSn-метод 6,4046 5,1633 4,7803 4,6092

DSnt-метод с DSn-распред 4,6704 4,6478 4,6450 4,6446

DSnt-метод с aSn-распред 3,6499 4,3770 4,5767 4,6523

Непрерывный Sn-метод [17] 4,4366 4,5859 4,6291 4,6467

Результаты расчетов трехобластной задачи, приведенные в таблице 3, показывают, что ОБ^-метод с ОБп-распределением сходится быстрее и превосходит по точности остальные методы, а ББтй-метод с аБп-распределением точнее классического ОБп-метода Полученные результаты позволяют заключить, что новые ОБ^-методы не уступают по точности ОБп-методу, они могут быть использованы для решения различных задач на определение критических параметров

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

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

Заключение

Результаты и выводы диссертационной работы

• ОБп-метод преобразован в две последовательно выполняемые схемы Первая - это разностная схема для ОДУ относительно полного потока на неосвещенных гранях, вторая - это схема его распределения по двум неосвещенным граням В новой двухэтапной форме классический ОЗп-метод представляет собой схему расщепления по причинам, вызывающим его теоретическую и практическую немонотонность

• Показано, что у ОБп-метода в новой двухэтапной форме схема для ОДУ аналогична простейшей разностной схеме ХТ-метода, а схема распределения эквивалентна дополнительным ап-проксимационным соотношениям ОБо-метода

• Сформулирован новый подход к построению разностных схем на Бп-сетках Суть подхода заключается в том, что разностная схема первого этапа заменяется в каждой ячейке на ОДУ баланса относительно полного потока с непрерывным изменением аргумента в самом ОДУ, в независимом источнике и источнике вторичных нейтронов В работе такое ОДУ построено

• Обоснована необходимость применения монотонных или квазимонотонных разностных схем 2-го порядка точности для численного решения ОДУ в каждой отдельной Бп-ячейке В этом случае первая причина немонотонности ОБп-метода, обусловленная аппроксимацией члена поглощения частиц стИ, полностью устраняется

• Предложена экономичная монотонная 2-го порядка точности разностная схема для решения ОДУ в Бп-ячейке и алгоритмы распределения полного потока по неосвещенным граням, что в совокупности образует новый численный метод - ОБп-метод характеристических трубок -ОБт-метод Эффективность новых ОБМ-схем подтверждена численными исследованиями различных классов задач переноса и кинетики нейтронов

Цитируемая в реферате литература

1 Б Карлсон, Дж Белл Решение транспортного уравнения Бп-методом - В сб «Физика ядерных реакторов». М, Атомиздат, 1959, стр 408—432

2 В Я Гольдин Методы расчета переноса нейтронов и горения в термоядерном изделии (1948-1960 гг) - Международный симпозиум, Дубна, 14-17 мая 1996 г В сб «Наука и общество история советского атомного проекта (40—50-е годы)», 1999, т 2, стр 497—501

3 Рихтмайер Р, Мортон К Разностные методы решения краевых задач M «Мир», 1972,

418 с

4 В С Владимиров Численное решение кинетического уравнения для сферы Вычислительная математика, 3, 1958, 3—33

5 Карлсон Б Г, Латроп К Д Теория переноса Метод дискретных ординат В сб Вычислительные методы в физике реакторов Под ред X Гринпсена, К Келбера и Д Окрента M, Атомиздат, 1972, стр 102—157

6 А В Никифорова, В А Тарасов, В Е Трощиев О решении кинетических уравнений дивергентным методом характеристик-ЖВМ и МФ, 1972,12, N 4, с 1041—1048

7 В Е Трощиев О математических свойствах Sn-методов решения кинетических уравнений-ЖВМ и МФ, 1975, 15, N5, с 1209—1221

8 В Е Трощиев Решение кинетического уравнения и уравнений квазидиффузии по согласованным разностным схемам - Численные методы решения задач математической физики (дополнение к ЖВМ и МФ, 6, № 4) M Наука, 1966 С 177—185

9 Л П Басс, A M Волощенко, Т А Гермогенова Методы дискретных ординат в задачах о переносе излучения ИПМ АН СССР, M, 1986

10 Carlson BGA method of characteristics and other improvements in solution methods for the transport equation Nuclear science and engineering 61,408—425(1976)

11 В А Елесин, В E Трощиев, В Ф Юдинцев, В И Федянин - Численная методика и организация программы для решения многогруппового нестационарного кинетического уравнения Сб Комплексы программ математической физики, Новосибирск, ВЦ СО АН СССР, 1972, 18—23

12 В Я Гольдин, А В Колпаков, А В Мисюрев Решение нестационарного уравнения переноса без явного выделения фронта - Препринт ИПМ АН СССР N 68, 1983

13 А Д Гаджиев, И А Кондаков, В H Писарев, О И Стародумов, А А Шестаков Метод дискретных ординат с искусственной диссипацией для численного решения уравнения переноса нейтронов // Вопросы атомной науки и техники Сер Математическое моделирование физических процессов 2003 Вып 4 с 13—24

14 О С Широковская Об одной модификации DSn-метода // Вопросы атомной науки и техники Сер Математическое моделирование физических процессов 1989 Вып 1 с 24—29

15 С В Мжачих, Е В Трошев, В Ф Юдинцев О некоторых свойствах DsJ¡ -схем для сферически-симметричного уравнения переноса // Вопросы атомной науки и техники Сер Математическое моделирование физических процессов 2000 Вып 2 с 21—31

16 В С Владимиров Математические задачи односкоростной теории переноса частиц -Тр Матем ин-таАНСССР, 1961,-158 с

17 О В Бутнева, Ю M Матвеев Численное сравнение дискретной и непрерывной аппроксимаций Sn-метода решения сферически-симметричного уравнения переноса ВАНТ, Серия Методики и программы численного решения задач математической физики, 1983, вып 3, С 36—40

18 В Е Трощиев, А В Ннфанова, Ю В Трощиев Характеристический подход к аппроксимации законов сохранения в кинетических уравнениях переноса излучений ДАН 2004, т 394, N 4, стр 454—458

19 Е В Грошев, А М Пастушенко, В Ф Юдинцев Об одной трехточечной разностной схеме с весовым множителем для уравнения переноса - Вопросы атомной науки и техники Сер Методики и программы численного решения задач математической физики 1985 Вып 2 с 87—96

20 W Н Reed New Difference Schemes for the Neutron Transport Equation Nucl Sei Eng, 46, 1971, p 309—315

21 С P Меркулова, В E Трощиев Монотонные разностные схемы для уравнения переноса и метод их построения Препринт ИАЭ им И В Курчатова, N5458/16, М, 1992

22 В Е Трощиев, А В Нифанова Построение и исследование разностных схем для уравнения переноса первого и второго порядка в плоском слое Препринт ТРИНИТИ N0052-A,(1999), 6 с

23 В Е Трощиев, Ю В Трощиев Монотонные разностные схемы с весом для уравнения переноса в плоском слое - Математическое моделирование, т 15, № 1,2003, с 3—13

Основные публикации по теме диссертации

1 Трощиев В Е, Нифанова А. В , Трощиев Ю В Характеристический подход к аппроксимации законов сохранения в кинетических уравнениях переноса излучений ДАН 2004, т 394, N 4, стр 454-458

2 В Е Трощиев, А В Нифанова Подход характеристических трубок к анализу DSn-метода и построение новых разностных схем на Sn-сетках - Математическое моделирование, т 18, № 7,2006, с 24—42

3 В Е Трощиев, А В Нифанова Построение и исследование разностных схем для уравнения переноса первого и второго порядка в плоском слое Препринт ТРИНИТИ N0052-A, (1999), 6 с

4 А В Нифанова, В Е Трощиев Нелинейная монотонная схема типа DSn-метода для уравнения переноса Научная сессия МИФИ-2001, сборник научных трудов, т 7, М, 2001, с 85—85

5 А В Нифанова, В Е Трощиев Обобщение метода характеристических трубок на Sn-сетки для сферически-симметричного уравнения переноса. Препринт ТРИНИТИ N0097-A, (2002), 16 с

6 А В Нифанова Характеристический подход к аппроксимации сферического уравнения переноса Научная сессия МИФИ-2003, сборник научных трудов, т 7, М, 2003, с 107—108

7 А В Нифанова, В Е Трощиев О методах распределения полного потока в схемах для сферического уравнения переноса Научная сессия МИФИ-2004, сборник научных трудов, т 7, М , 2004, с 98—99

8 В Е Трощиев, А В Нифанова Характеристический анализ непрерывного Sn-метода Научная сессия МИФИ-2007, сборник научных трудов, т 7, М , 2007, с 100—101

НПФАИОВА Александра Васильевна

^РАКТЕРИСТИЧЕСКИЕ Бп-МЕТОДЫ ДЛЯ КИНЕТИЧЕСКОГО УРАВНЕНИЯ ПЕРЕНОСА НЕЙТРОНОВ

В СФЕРИЧЕСКИХ СИСТЕМАХ

Автореферат диссертации

Подписано в печать //. 0?. 2008 Формат 60x90/16 Гарнитура Times New Roman Уч изд. листов 1, тираж 70 экз

Заказ №217

Подготовлено к изданию и отпечатано в ФГУП «ГНЦ РФ ТРИНИТИ», 142190, г Троицк Моек обл.ул Пушковых, влад 12

Введение.

Глава 1. Кинетическое уравнение переноса нейтронов для сферы и математические постановки задач о его решении.

§1 Интегро-дифференциальное уравнение переноса в частных производных первого порядка.

§2 Семейство характеристик уравнения переноса в фазовом пространстве (г, //) и семейство обыкновенных интегро-дифференциальных уравнений (ОИДУ) на характеристиках.

§3 Классы задач для уравнений переноса pi итерационные процессы их решения.

Комментарии к главе

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

§1 Разностная схема DSn-метода в ячейке баланса.

§2 Преобразование разностной схемы DSn-метода в схему для ОДУ относительно полного потока и схему распределения потока по неосвещенным граням Sn-ячейки (расщепление DSn-схемы).

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

§4 Дискретная форма и расщепление непрерывного Sn-метода на схему для ОДУ и схему распределения полного потока.

§5. О расщеплении балансных разностных схем для нестационарных и многомерных кинетических уравнений переноса частиц.

Выводы и комментарии по главе 2.

Глава 3. Обобщение метода характеристических трубок на ячейки-трубки общего вида и построение двухэтапных DSn-методов квазимонотонного типа.

§1 Построение осредненного ОДУ баланса в характеристических ячейках-трубках общего вида.

§2 Sn-ячейка как ячейка-трубка. $3 Положительные монотонные и квазимонотонные схемы 2-го порядка точности для ОДУ *.

§4 Схемы распределения полного потока.

Выводы к главе 3.

Глава 4. Численное исследование точности двухэтапных характеристических DSn-методов для различных классов сферических задач.

51 Тестирование новых схем на решении задач с источником.

52 Применение двухэтапных методов к решению задач на определение критических параметров.

Выводы к главе 4.

Кинетические уравнения переноса излучений - нейтронов, фотонов и других частиц [1-3] являются математическим описанием физических законов сохранения. Кратко, без детализации, эти законы можно свести к двум положениям - это сохранение числа частиц при их движении без взаимодействия со средой (дивергентность уравнения переноса) и определенное количественное соответствие между числом исчезнувших частиц при их взаимодействии со средой и числом вновь появившихся вторичных частиц (самосогласованность столкновительных членов в кинетическом уравнении). Точное выполнение законов сохранения в дискретных моделях переноса и кинетики частиц {консервативность модели) является одним из наиболее важных аппроксимационных требований. Поэтому при построении консервативных сеточных аппроксимаций методом баланса [4] кинетические уравнения переноса частиц [1-3] принято рассматривать в дивергентной форме ([5-18] и др.):

LN(r,&) = div(N • п)+ a(r)N = £>{?,&) qI?,Q)= ТЫ(гД)+ G(r,n), TN(r,&)= |{/(й' -> а)\ф,й')ь/а' (EU) и строго выполнять сеточную дивергентность и сеточную самосогласованность столкновительных членов cr(r)N, TN и таким образом обеспечивать консервативность модели в дискретной форме.

Следующее важное требование, предъявляемое к сеточным аппроксимациям кинетических уравнений (В. 1), заключается в построении и использовании разностных схем, порядок точности которых не ниже второго на гладких решениях (это требование обусловлено очень низкой точностью схем первого порядка). Так как уравнение переноса (В.1) линейное, то невозможно построить линейную разностную схему второго порядка точности, которая была бы теоретически положительной и монотонной [19]. Немонотонность схем второго порядка точности часто проявляется и в практических расчетах, поэтому вопрос монотонизации сеточных решений для этого класса схем с обеспечением хорошей точности расчетов является исключительно актуальным. Частичным решением проблемы немонотонности сеточных решений для схем второго порядка точности является разработка и локальное по ячейкам сетки применение «алгоритмов коррекции потока» (АКП) [12]. Однако в АКП очень трудным вопросом является математически строгое определение ячеек сетки, в которых проявляется немонотонность, и оценка ошибки в сеточном решении. Достаточно простого и общего решения этого вопроса не существует.

Еще одно важное требование к численным методам решения уравнений переноса заключается в необходимости построения схем бегущего счета (СБС), то есть схем с треугольными матрицами у полной системы сеточных уравнений. Это требование связано с тем, что уравнение (В.1) является интегральным относительно процесса рассеяния вторичных частиц по направлениям полета, что в свою очередь требует применения итерационных алгоритмов по интегралу TN{r,£lj в кинетическом уравнении (В.1). По своему смыслу интеграл TN собирает на направление Q частицы, рассеянные со всех направлений О! на единичной сфере с центром в точке 7. Поэтому реальным способом численного решения задач переноса и кинетики нейтронов является метод итераций по интегралу рассеяния вторичных нейтронов

V v-1 v-1 V-1

LN = Q, Q=TN+G. (В.2)

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

LhNh = Qh, v = 1,2,3,. (В.З) необходимо строить разностные схемы с треугольными матрицами у полной системы сеточных уравнений (схемы бегущего счета). В работе [20] показано, что в общем случае такие схемы можно строить на сетках, состоящих из выпуклых ячеек соответствующей размерности. Схемами бегущего счета, в частности, являются Sn-схемы [21-23,7,10,11], схемы метода характеристик [24,11] и характеристических трубок (ХТ-схемы) [9], многомерные схемы [13,15,16,18].

Для широкого круга задач нейтронно-ядерной физики со сферической симметрией процессов переноса и кинетики нейтронов основными математическими методами их численного решения являются Sn-методы и методы характеристик. Первые варианты этих методов были независимо сформулированы в конце 1940-х, начале 1950-х годов в работах по атомным проектам США и СССР. Это Sn-метод Карлсона [22] и КН-схема Гольдина [6], метод прямого интегрирования Рихтмайера [25] и метод характеристик Владимирова [24]. В 1960-х годах были предложены дискретный Sn-метод (DSn-метод) [22] и метод характеристических трубок (ХТ-метод) [9], которые представляют собой развитие и обобщение в определенных направлениях первоначальных методов. Эти два метода и их модификации наиболее широко используются в практических расчетах задач нейтронно-ядерной физики и физики взаимодействия теплового излучения с веществом. К таким задачам относятся определение функции плотности нейтронов и интенсивности излучений в стационарных и нестационарных задачах, расчет критических параметров в сферических системах, расчеты потоков в задачах защиты от излучений и другие. Достоинства этих методов и имеющиеся в них недостатки обуславливают актуальность их дальнейшего развития и совершенствования.

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

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

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

Ячейками сетки в DSn-методе являются прямоугольники в фазовом пространстве координат г (0<r<R) и направлений полета нейтронов (J. (-1<ц<1) (рис.2.1.1, рис. 2.1.2). Уравнение баланса нейтронов в каждой ячейке сетки представляет собой аппроксимацию уравнения в частных производных 1-го порядка по переменным ги ц. Ячейками сстки в XT-методе являются отрезки характеристических трубок - это криволинейные четырехугольники, образованные характеристиками-траекториями в фазовом пространстве (r,|i) и линиями сетки по переменной г (рис. 2.3.2). При этом в каждой ячейке-трубке выполняется осредненное ОДУ баланса нейтронов.

DSn- и XT-методы - принципиально разные по своей сути, но их первоначальной аппроксимационной основой являются сеточные уравнения баланса, записанные соответственно на Sn- и Т-сетках. Это аппроксимационное свойство в определенной степени их сближает.

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

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

Sn-метод был сформулирован Б. Карлсоном для сферически-симметричного уравнения переноса в недивергентной форме, в связи с чем для выполнения сеточной дивергентности Sn-схемы потребовались непростые дополнительные преобразования расчетных формул [21]. В Sn-методеиспользуется прямоугольная сетка по переменным г и М (о < г < R,-l < pi < /)> которую мы будем называть Sn-сеткой. Сеточная функция в пределах отдельной прямоугольной ячейки предполагается линейной по каждой переменной (непрерывный Sn-метод). В последующих работах Б. Карлсоном был сформулирован дискретный Sn-метод (DSn-метод) [22], в котором в качестве сеточной функции используются средние значения потоков частиц на гранях (сторонах) Sn-ячейки, а уравнение баланса замыкается еще одним уравнением, связывающим значения сеточной функции на 4-х гранях Sn-ячейки.

В.Я. Гольдиным при построении КН-схемы (консервативная нестационарная схема) была использована дивергентная форма уравнения переноса, что автоматически обеспечивало ее консервативность [5,6]. Позже была показана эквивалентность КН-схемы и схемы Sn-метода для уравнения в дивергентной и недивергентной формах, получены оценки порядка аппроксимации схемы в окрестности и вне центра сферы [7, 8, 10]. К настоящему времени установлены многие математические свойства Sn-методов, из которых отметим следующие.

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

2. Sn-методы — это линейные схемы второго порядка аппроксимации на гладких решениях, которые, согласно теореме Годунова, являются теоретически (и практически) немонотонными.

3. В окрестности центра сферы у Sn-метод а на единицу понижается порядок аппроксимации и точности, а DSn-метод имеет нулевой порядок аппроксимации и точности [10].

4. В окрестности центра сферы при измельчении сетки по радиусу решение сходится к const по угловой переменной без ее измельчения [10].

5. Для уравнений квазидиффузии существует разностная схема, согласованная с DSn-методом [7].

В дальнейшем, в ряде работ [13-18] и др. были предложены конечно-разностные схемы для многомерных уравнений переноса в криволинейных системах координат, которые можно рассматривать как обобщения DSn-метода на многомерные косоугольные сетки.

В работе [26] было рассмотрено семейство разностных схем типа Sn-метода, которые различаются аппроксимацией (непрерывной или дискретной) потоковых и столкновительных членов. Сравнение численных результатов, проведенное авторами, показало, что практическая точность непрерывного Sn-метод с непрерывной же аппроксимацией столкновительных членов выше, чем точность DSn-метода, но он требует примерно в 2,5 раза большего числа арифметических операций. Непрерывный Sn-метод иногда используется для построения схем типа DSn-метода с улучшенными свойствами [12].

Существуют подходы (они берут начало из работ [ 10,27]), в которых предлагается расширенная трактовка DSn-метода, а именно, он рассматривается как представитель некоторого, например, четырехпараметрического, как в работах [28], семейства схем. В работе [28а] такое семейство схем строится на основе предложенного авторами аппарата взаимных преобразований дивергентной и недивергентной разностных форм оператора переноса, а также рассматриваются алгоритмы подбора параметров схемы, и некоторые из этих алгоритмов используют соотношения непрерывного Sn-метода. Эффективность метода авторы подтверждают расчетами сферически-симметричной задачи (однородный шар) с параметрами близкими к критическим. Предложенные авторами схемы дают заметное по сравнению с DSn-методом повышение точности на грубых сетках.

Немонотонность является существенным недостатком DSn-метода (равно как и Sn-метода). Уравнение переноса может быть многомерным даже для пространственно одномерных задач, поэтому расчеты нередко ведутся на грубых сетках. Низкая точность схем первого порядка не позволяет их применять для получения монотонных решений с удовлетворительной точностью. Требуются схемы второго или более высокого порядка точности. В обзоре [12] подчеркивается, что наиболее сложно обеспечить одновременное требование положительности и высокой точности. Поэтому практически одновременно с Sn-методом начали появляться методы его коррекции. В качестве коррекции Sn-схем одним из первых был предложен метод балансного зануления отрицательных потоков в ячейке [29], который, несмотря на имеющиеся недостатки, успешно используется в некоторых задачах. Использовался также пересчет ячейки, в которой возникли отрицательные потоки, по положительной схеме первого порядка (шаговая или St-схема). Как отмечается в обзоре [12], такие методы коррекции могут ухудшить сходимость итерационного процесса в сложных задачах. Более успешно применяется взвешенная схема. Она отличается от Sn-схемы 2-го порядка точности тем, что в дополнительные соотношения схемы вводится весовой множитель. Взвешенная схема занимает некоторое промежуточное положение между схемами 1-го и 2-го порядка аппроксимации.

Для построения монотонных схем для разных классов уравнений неоднократно привлекался характеристический подход [30-32] и др. В работе [30] была предложена нелинейная схема с монотонной аппроксимацией дифференциального оператора для модельного уравнения переноса. Искомое значение в неосвещенной вершине прямоугольной ячейки вычислялось в результате параболической интерполяции со сглаживанием. Для построения интерполяции брались значения плотности потока в точках на одном из ребер ячейки и его продолжении. В эти точки попадали характеристики дифференциального оператора, выпущенные из освещенных вершин, а значения плотности потока в них определялись точным решением уравнения переноса вдоль выпущенных характеристик. Сама задача при этом решалась методом квазидиффузии [33]. Этот подход получил развитие в работе [34]. В двумерном уравнении интерполяция строилась на произвольной прямой, проведенной через неосвещенную вершину ячейки, а в многомерном случае на произвольной плоскости, проходящей через неосвещенную вершину ячейки.

В работе [35] был предложен метод «мягкой коррекции», в котором весовая функция в конкретной ячейке записывается на основе характеристического представления уравнения переноса в этой ячейке и зависит от ее оптической толщины. Этот метод коррекции дал начало адаптивным схемам [12, 36], в которых учитывается структура решения в областях с различными физическими параметрами. Некоторые из этих методов, особенно последние, достаточно эффективны при определенных условиях, но не являются достаточно общим решением проблемы немонотонности, а их обобщения на многомерные уравнения очень сложны.

В работе [37] отмечается, что нефизическое осциллирующее решение в оптически плотной среде DSn-метод дает из-за линейной интерполяции столкновительного члена, описывающего поглощение частиц в кинетическом уравнении. В связи с этим авторами предлагается DSn-метод с добавкой диссипативных членов, который они называют DDAD-схемой (Diamond Difference Artificial Dissipation). В работе приводятся расчеты гетерогенных слоистых задач (для плоской, сферически-симметричной и двумерной цилиндрической геометрии), которые были выполнены по ромбовой (алмазной) схеме (DD) с пересчетом отрицательных значений по шаговой схеме (DD/St) и по DDAD-схемс с пересчетом отрицательных значений по шаговой схеме (DDAD/St). Результаты расчетов по новой схеме заметно лучше, чем по DD/St-схеме. Этот же алгоритм DDAD применяется и в работе [36] для численного решения двумерного уравнения переноса теплового излучения в многогрупповом приближении.

Как отмечается в работе [14], в DSn-методе аппроксимация уравнения переноса по угловой переменной даже в случае изотропного рассеяния не всегда удовлетворительна, особенно в задачах о прохождении частиц на большие расстояния. Автор предлагает нелинейную модификацию дополнительных соотношений и квадратурных формул DSn-метода, что позволяет получать решение с хорошей точностью на довольно грубой сетке в пространстве направлений. Эта же проблема подробно рассматривается в работе [39]. Чтобы выделить ошибки, порождаемые разностной аппроксимацией по угловой переменной автор вместо разностной аппроксимации по пространственной переменной (радиусу в сферически-симметричной задаче) предлагает решение системы связанных друг с другом обыкновенных дифференциальных уравнений с использованием стандартного пакета.

В ряде работ [12,36,40-43], были предложены схемы более высокого порядка аппроксимации, чем DSn-метод, некоторые из них обладают и лучшими свойствами в плане положительности и монотонности, но при этом они сложнее для численной реализации, и особенно для обобщения на многомерные геометрии. Задача построения монотонных положительных разностных схем для уравнений переноса с аппроксимацией не ниже второго порядка, достаточно простых и универсальных, остается актуальной.

В 1951-1952 г.г. B.C. Владимировым был разработан метод характеристик для численного решения задач на собственные значения и задач с источником в многослойных сферически-симметричных системах [24,44-46,11]. Идея метода состояла в том, что в результате замены переменных интегро-дифференциальное уравнение с частными производными сводится к уравнениям вдоль характеристик, то есть вдоль траекторий полета частиц. Решение можно находить, если ввести сетку * на характеристиках, то есть разбить их на интервалы, внутри каждого интервала приблизить правую часть уравнения линейной функцией и произвести точное интегрирование последовательно на всех интервалах. Такой способ всегда дает положительное решение. Недостатком метода является его неконсервативпость, которая приводит во многих практически важных задачах к снижению точности численного решения.

Метод характеристических трубок — ХТ-метод [9,12], предложенный В.Е. Трощиевым (кандидатская диссертация, 1966г.), является продолжением и развитием метода характеристик Владимирова. В сферически-симметричной геометрии характеристическая трубка - это область, заключенная между двумя «соседними» характеристиками дифференциального оператора уравнения переноса, которая разбивается на Т-ячейки баланса прямыми постоянных радиусов (торцы Т-ячейки). В каждой Т-ячейке вводится сеточная функция на торцах и записывается разностное уравнение баланса частиц, которое преобразуется таким образом, что уже может рассматриваться как разностная аппроксимация обыкновенного дифференциального уравнения. Важным свойством ОДУ является то, что его можно рассматривать как относительно функции плотности потока частиц вдоль средней характеристики, так и функции потока вдоль трубки. В последнем случае это есть ОДУ баланса частиц в Т-ячейке. Для решения ОДУ существует ряд методов, дающих положительное и монотонное решение. Таким образом, XT-метод положителен, монотонен и имеет 2-й порядок точности, в нем точно выполняются законы сохранения нейтронов. Широкое применение XT-метода для решения задач о переносе нейтронов и у-квантов [9] показало его хорошие практические качества и подтвердило его теоретические свойства. Недостатком XT-метода является сложность его расчетной сетки и принципиальные трудности его обобщения на нестационарные и многомерные задачи ( в отличие от DSn-метода).

В работе [46] был предложен новый подход к построению разностных схем для многомерных уравнений переноса, который является нетривиальным обобщением ХТ-метода на выпуклые многомерные косоугольные сетки. Суть подхода [46] состоит в трактовке любого типа сеточных ячеек как характеристических трубок с двумя многогранными торцами - освещенным и неосвещенным - и последующем построении осредненного ОДУ баланса частиц в такой ячейке-трубке. В итерационном процессе по интегралу столкновений ОДУ относительно полного потока нейтронов вдоль ячейки-трубки может быть эффективно решено по консервативным положительным, монотонным или квазимонотонным нелинейного типа схемам 2—го порядка точности. Такие схемы построены в работах [47-51, 12], причем схемы [51] являются линейными. После решения ОДУ возникает задача консервативного распределения найденного полного потока по неосвещенным граням исходной ячейки. Задача распределения уже известного полного потока становится самостоятельной и должна решаться для каждой конкретной схемы с учетом структуры сеток и способов задания сеточных функций. В частности, распределение может быть выполнено на основе дополнительных интерполяционных соотношений, применяемых в DSn-методе и его обобщениях на многомерные уравнения и сетки.

В данной работе характеристический подход [46] развит для сферически-симметричного уравнения переноса частиц на сетках Sn-методов [52-55]. В этом случае особенности формулировок и исследования характеристического подхода обусловлены криволинейностью сферической системы координат, что наглядно проявляется при анализе дискретного Sn-метода.

Целью работы было построение новых разностных схем на базе схем DSn- и ХТ-методов на Sn-сетках (DSn-методы характеристических трубок или DSnt-методы) и численное исследование этих схем на монотонность и точность в расчетах основных классов задач переноса и кинетики нейтронов.

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

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

Во второй главе рассматриваются Sn- и XT-методы решения задач переноса. Дискретный и непрерывный Sn-методы подвергаются анализу с точки зрения характеристического подхода [46, 55, 56]. Введением дополнительной сеточной функции — полного потока частиц — DSn-метод преобразуется в двухэтапную схему - схему метода характеристических трубок для полного потока и DSn-схему его распределения по неосвещенным граням Sn-ячейки, и таким образом, расщепляется по причинам, вызывающим его немонотонность. В главе также обсуждается вопрос о возможности расщепления балансных разностных схем для многомерных кинетических уравнений переноса частиц.

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

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

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

Основные результаты работы следующие:

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

2. Математические понятия инварианта переноса и среднего расстояния, ранее введенные в методе характеристических трубок, обобщены на сетки произвольной формы. Установлена их связь с фазовым объёмом сеточных ячеек в сфере. В Sn-ячейке, трактуемой как характеристическая трубка, построено ОДУ баланса относительно полного потока и функция независимого источника с непрерывным изменением аргумента - расстояния от освещенных граней до неосвещённых.

3. Для ОДУ в Sn-ячейке-трубке предложена экономичная монотонная 2-го порядка точности разностная схема для определения полного потока. Эта схема и различные алгоритмы распределения полного потока по неосвещенным граням представляют собой новые численные методы (DSn-методы характеристических трубок - DSnt-методы), в которых полностью устранена причина немонотонности, обусловленная DSn-аппроксимацией столкновительных членов в уравнении переноса и кинетики нейтронов.

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

Результаты работы докладывались и обсуждались на Научных сессиях МИФИ (2001, 2003, 2004, 2007), на семинаре «Нейтроника-2005» в Обнинске, на семинаре В. Я. Гольдина Института математического моделирования РАН, на семинаре В.М.Головизнина ИБРАЭ РАН.

По теме диссертации опубликовано 8 печатных работ, из них две статьи в реферируемых журналах - «Доклады академии наук» и «Математическое моделирование», два препринта ГНЦ РФ ТРИНИТИ и четыре сообщения в сборниках трудов Научных сессий МИФИ.

ВЫВОДЫ

• DSn-метод преобразован в две последовательно выполняемые схемы. Первая — это разностная схема для ОДУ относительно полного потока на неосвещенных гранях, вторая — это схема его распределения по двум неосвещенным граням. В новой двухэтапной форме классический DSn-метод представляет собой схему расщепления по двум причинам, обусловливающим его теоретическую и . практическую немонотонность.

• Показано, что у DSn-метода в новой двухэтапной форме схема для ОДУ аналогична простейшей разностной схеме метода характеристических трубок — XT-метода, а схема распределения эквивалентна дополнительным интерполяционным соотношениям дискретного Sn-метода - DSn-метода.

• Сформулирован новый подход к построению разностных схем на Sn-сетках. Суть подхода заключается в том, что разностная схема первого этапа заменяется в каждой ячейке на ОДУ баланса относительно полного потока с непрерывным изменением аргумента в самом ОДУ, в независимом источнике и источнике вторичных нейтронов. В работе такое ОДУ построено.

• Применение положительных монотонных или квазимонотонных разностных схем 2-го порядка точности для численного решения ОДУ в каждой отдельной Sn-ячейке полностью устраняет первую причину немонотонности DSn-метода, обусловленную аппроксимацией члена поглощения частиц в кинетическом уравнении.

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

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

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

ТРИНИТИ Ежову А. А. за помощь в работе.

ЗАКЛЮЧЕНИЕ

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

Выполняемая первой схема — это схема метода характеристических трубок (ХТ-метод [9]) для некоторого осредненного ОДУ относительно полного потока частиц. По этой схеме в DSn-методе рассчитываются физические процессы поглощения нейтронов при их взаимодействии со средой и полный поток нейтронов через неосвещенные грани Sn-ячеек. Таким образом, схема первого этапа — это схема, определяющая баланс нейтронов в DSn-методе.

Выполняемая второй схема - это схема распределения найденного полного потока по двум неосвещенным граням. Эта схема в точности эквивалентна дополнительному интерполяционному соотношению DSn-метода, которое вместе с сеточным уравнением баланса DSn-метода дает замкнутую систему двух уравнений для искомых значений сеточной функции в рассчитываемой Sn-ячейке. Схема DSn-распределения полного потока — это аппроксимация и расчет процесса переноса нейтронов в Sn-ячейке без их взаимодействия со средой по известным значениям полных потоков. Таким образом, схема второго этапа - это схема решения кинетического уравнения в Sn-ячейке, трактуемой как полость.

Расщепление DSn-метода на две схемы — это новое свойство метода, которое является наиболее важным на пути построения новых схем — DSn-схем характеристических трубок (DSnt-схемы). Свойство «расщепления» также показывает глубокую взаимную связь Sn-методов и методов характеристик, которые рассматривались ранее как совершенно разные методы, можно сказать, как альтернативные. Расщепление означает отделение друг от друга двух причин, вызывающих немонотонность DSn-метода. Первая причина — это линейная 2-го порядка точности аппроксимация процессов кинетики частиц в столкновительных членах кинетического уравнения. Вторая причина — это аппроксимация процесса бесстолкновительного переноса частиц в уравнении первого порядка с частными производными.

На основе проведенного анализа DSn-метода (расщепления) в работе развит характеристический подход (подход характеристических трубок) к построению консервативно-характеристических схем на прямоугольных Sn-сетках для сферически-симметричного кинетического уравнения переноса. Главная суть XT подхода — это построение в Sn-ячейках ОДУ баланса частиц относительно функции полного потока с непрерывным изменением аргумента — расстояния от неосвещенной поверхности Sn-ячейки до освещенной и решение этого ОДУ по монотонным схемам второго порядка точности. Это первый этап в построении и реализации новых разностных схем. На этом этапе в новых схемах полностью устраняется первая причина немонотонности Sn-методов.

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

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

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

1. Дэвисон Б. Теория переноса нейтронов. М., Атомиздат, 1960.

2. Я.Б. Зельдович, Ю.П. Райзер. Физика ударных волн и высокотемпературных гидродинамических явлений. М. Наука, 1966, 686с.

3. Марчук Г. И., Лебедев В. И. Численные методы в теории переноса нейтронов. М., Атомиздат, 1981,454 с.

4. А.А. Самарский. Теория разностных схем. М., Наука, 1983

5. Самарский А.А. Прямой расчет мощности взрыва. Международный симпозиум, Дубна, 14-17 мая 1996 г. В сб.: "Наука и общество: история советского атомного проекта (40-е -50-е годы)", 1997, том. 1, стр. 214-222.

6. В.Е. Трощиев. Решение кинетического уравнения и уравнений квазидиффузии по согласованным разностным схемам. Численные методы решения задач математической физики (дополнение к ЖВМ и МФ, 6,№ 4). М.: Наука, 1966. С.177-185.

7. В.Е. Трощиев, В.Ф. Юдинцев, В.И. Федянин. Об ускорении сходимости итераций при решении кинетического уравнения. ЖВМ и МФ, 1968, 8, №2, с.452-458.

8. А.В. Никифорова, В.А. Тарасов, В.Е. Трощиев. О решении кинетических уравнений дивергентным методом характеристик.-ЖВМ и МФ, 1972, 12, N4,с. 1041-1048.

9. В.Е. Трощиев. О математических свойствах Sn-методов решения кинетических уравнений.-ЖВМ и МФ, 1975, 15, N5, с.1209-1221.

10. Смелов В.В. Лекции по теории переноса нейтронов. М., Атомиздат, 1978.

11. Басс, A.M. Волощенко, Т.А. Гермогенова. Методы дискретных ординат в задачах о переносе излучения. ИПМ АН СССР, М., 1986.

12. В.Е. Трощиев, В.А. Шумилин. Разностная схема решения двумерного уравнения переноса на нерегулярных четырехугольных сетках. ЖВМ и МФ, 1986, 26,N2, с.230-241.

13. О. С. Широковская. Об одной модификации DSn-метода. // Вопросы атомной науки и техники. Сер. Математическое моделирование физических процессов. 1989. Вып.1. с.24-29.

14. Л.П. Федотова, P.M. Шагалиев. Конечноразностный КМ-метод для двумерных нестационарных процессов переноса в многогрупповом кинетическом приближении. Математическое моделирование. 1991. т.3,№6. С.29-41.

15. В.А. Елесин, В.Е. Трощиев, В.Ф. Юдинцев. Развитие численных методов и программ расчета одномерных спектральных задач переноса теплового излучения во ВНИИЭФ. — ВАНТ, Серия: Математическое моделирование физических процессов, 2002, Вып. 1, с. 11-28.

16. А. Н. Москвин, В. А. Шумилин. Методика решения двумерного уравнения переноса на нерегулярных многоугольных сетках. // Вопросы атомной науки и техники. Сер. Математическое моделирование физических процессов. 2005. Вып.1. с.31-40.

17. Калиткин Н. Н. Численные методы. М., Наука, 1978.

18. В.Е. Трощиев. О классах сеток, допускающих консервативные аппроксимации двумерного оператора переноса треугольным разностным оператором. ЖВМ и МФ, 1976, 16,№ 3, с.793-797.

19. Б. Карлсон, Дж. Белл. Решение транспортного уравнения Sn-методом. В сб. "Физика ядерных реакторов". М., Атомиздат, 1959, стр.408—432.

20. Б. Карлсон. Численное решение задач кинетической теории нейтронов. В сб. "Теория ядерных реакторов". М., Госатомиздат, 1963, 243-258.

21. Карлсон Б.Г., Латроп К.Д. Теория переноса. Метод дискретных ординат. В сб.: Вычислительные методы в физике реакторов. Под ред. X. Гринпсена, К. Келбера и Д. Окрента. М., Атомиздат, 1972, стр.102-157.

22. B.C. Владимиров. Численное решение кинетического уравнения для сферы. Вычислительная математика, 3, 1958, 3-33.

23. Рихтмайер Р., Мортон К. Разностные методы решения краевых задач. М. «Мир», 1972, 418 с.

24. В.А. Елесин, В.Е. Трощиев, В.Ф. Юдинцев, В.И. Федянин Численная методика и организация программы для решения многогруппового нестационарного кинетическогоуравнения. Сб.: Комплексы программ математической физики, Новосибирск, ВЦ СО АН СССР, 1972, 18-23.

25. В. Я. Гольдин, Н. Н. Калиткин, Т. В. Шишова. Нелинейные разностные схемы для гиперболических уравнений. Журнал вычислительной математики и математической физики, 1965, т.5, №5, с.938-944.

26. Головизнин В. М. Балансно-характеристический метод численного решения уравнений газовой динамики. ДАН, 2005, т.403, №4, с. 1-6.

27. Аристова Е. Н., Гольдин В. Я., Дементьев А. С. Разностное решение двумерного стационарного уравнения переноса з переменных Владимирова. Математическое моделирование, 2006, т. 18, №6, с.44-52.

28. В. Я. Гольдин. Квазидиффузионный метод решения кинетического уравнения. ЖВМ и МФ, 1964, т.4, №6, с. 1078-1087.

29. В. Я. Гольдин, А. В. Колпаков, А. В. Мисюрев. Решение нестационарного уравнения переноса без явного выделения фронта. Препринт ИПМ АН СССР N68, 1983.

30. Carlson В. G. A method of characteristics and other improvements in solution methods for the transport equation. Nuclear science and engineering: 61, 408-425 (1976).

31. Т. А. Гермогенова, A. M. Волощенко. К развитию метода дискретных ординат. // Вопросы атомной науки и техники. Сер. Физика и техника ядерных реакторов. 1985. Вып. 5. С. 57.

32. К. D. Lathrop. A Comparison of Angular Difference Schemes for One-Dimensional Spherical Geometry SN Equations. Nuclear Science and Engineering: 134, 239-264 (2000).

33. E. В. Диянкова, О. С. Широковская. LD-схема для уравнений переноса в сферической геометрии. // Вопросы атомной науки и техники. Сер. Математическое моделирование физических процессов. 1989. Вып.1. с.40-43.

34. Е. В. Диянкова, О. С. Широковская. Разностная схема повышенного порядка аппроксимации для уравнения переноса. Математическое моделирование, 1994, т.6, №2, с.113-122.

35. Владимиров B.C. Математические задачи односкоростной теории переноса частиц. -Тр. Матем. ин-та АН СССР, 1961, 158 с.

36. Марчук Г.И. Численные методы расчета ядерных реакторов. М., Атомиздат,1958,381 с.

37. Трощиев В. Е., Нифанова А. В., Трощиев Ю. В. Характеристический подход к аппроксимации законов сохранения в кинетических уравнениях переноса излучений. ДАН, 2004, т394, N4, стр.454-458.

38. W. Н. Reed. New Difference Schemes for the Neutron Transport Equation. Nucl. Sci. Eng., 46, 1971, p.309-315

39. С. P. Меркулова, В. E. Трощиев. Монотонные разностные схемы для уравнения переноса и метод их построения. Препринт ИАЭ им. И.В. Курчатова, N5458/16, М., 1992.

40. В. Е. Трощиев, А. В. Нифанова. Построение и исследование разностных схем для уравнения переноса первого и второго порядка в плоском слое. Препринт ТРИНИТИ N0052-A,(1999), 6с.

41. А. В. Нифанова, В. Е. Трощиев. Нелинейная монотонная схема типа DSn-метода для уравнения переноса. Научная сессия МИФИ-2001 сборник научных трудов, т.7, М., 2001, с.85-85.

42. Трощиев В.Е., Трощиев Ю.В. Монотонные разностные схемы с весом для уравнения переноса в плоском слое. Математическое моделирование, т. 15, № 1, 2003, с. 3-13.

43. А.В. Нифанова. Характеристический подход к аппроксимации сферического уравнения переноса. Научная сессия МИФИ-2003 сборник научных трудов, т.7, М., 2003, с. 107-108.

44. А. В. Нифанова, В. Е. Трощиев. О методах распределения полного потока в схемах для сферического уравнения переноса. Научная сессия МИФИ-2004 сборник научных трудов, т.7, М., 2004, с.98-99.

45. А. В. Нифанова, В. Е. Трощиев. Обобщение метода характеристических трубок на Sn-сетки для сферически-симметричного уравнения переноса. Препринт ТРИНИТИ N0097-А,(2002), 16с.

46. В. Е. Трощиев, А. В. Нифанова. Подход характеристических трубок к анализу DSn-метода и построение новых разностных схем на Sn-сетках. Математическое моделирование, т. 18, № 7, 2006, с. 24-42.

47. В.Е. Трощиев, А.В. Нифанова. Характеристический анализ непрерывного Sn-метода. Научная сессия МИФИ-2007 сборник научных трудов, т.7, М., 2007, с. 100-101.

48. В. Я. Гольдин. Характеристическая разностная схема для нестационарного кинетического уравнения. ДАН СССР, 1960,т.133, №4, с.748-751.

49. В.Е. Трощиев. Метод построения блочно-треугольных разностных схем для уравнения переноса в самосопряженной форме. Математическое моделирование. 1998, 10, №1, с.117-125.