Содержание к диссертации
Введение
ГЛАВА 1. Модифицированная численная схема теплопроводного этапа комплекса программ НЗТ 21
1. Введение 21
2. Система дифференциальных уравнений модели 23
3. Разностная схема 25
4. Тестовый расчет 32
5. Параллельные аспекты схемы 39
6. Заключение 41
ГЛАВА 2. Алгоритм построения фронтов ударных волн и контактных границ с искусственной вязкостью 42
1. Расчет фронта ударной волны 43
1.1 Структура сетки в расчетной области 43
1.2 Алгоритм переинтерполяции граничной скорости из полуцелых точек в целые 45
2. Сравнение алгоритма построения фронтов ударных волн и кон тактных границ с искусственной вязкостью с другими алгоритмами 47
3. Тестовый расчет 1 49
3.1 Задача об обтекании сферы, постановка 49
3.2 Расчет по алгоритму с равномерной расстановкой узлов вдоль выделенной ударной волны 50
3.3 Расчеты по алгоритму с вязкостью 51
4. Рисунки 54
5. Апробация алгоритма построения фронтов ударных волн и контактных границ с искусственной вязкостью на задаче о вытекании вещества через торцы мишени 61
5.1 Постановка задачи 61
5.2 Результаты расчета 62
5.3 Рисунки 63
6. Выводы 65
ГЛАВА 3. Численные расчеты 66
1. Конструкция цилиндрической мишени с урановым пушером, сжи маемой безударным способом 67
1.1 Постановка задачи 68
1.2 Результаты расчета 69
1.3 Выводы 71
2. Сравнение двумерной методики НЗТ с одномерной методикой СНД 72
2.1 Постановка задачи 73
2.2 Результаты сравнения 73
2.3 Выводы 77
2.4 Рисунки 78
3. Методика Н2Т-1Р 87
3.1 Система уравнений методики НЗТ 88
3.2 ID система уравнений модели 2Т-1Р 91
3.3 Уравнение переноса 94
3.4 Алгоритм задачи 99
3.5 Постановка расчетной задачи 102
3.6 Параллельные аспекты 103
3.7 Расчеты 105
3.8 Выводы 111
ЗАКЛЮЧЕНИЕ 112
Литература
- Система дифференциальных уравнений модели
- Алгоритм переинтерполяции граничной скорости из полуцелых точек в целые
- Апробация алгоритма построения фронтов ударных волн и контактных границ с искусственной вязкостью на задаче о вытекании вещества через торцы мишени
- Сравнение двумерной методики НЗТ с одномерной методикой СНД
Введение к работе
Бурное развитие науки, в частности физики, в начале двадцатого века привело к появлению такого понятия, как ядерная энергия. Если учесть, что все предыдущие виды энергии были найдены обычными людьми и исследователями, то здесь за дело взялись ученые-физики. Уже в пятидесятые годы прошлого века ученые нашли способ управлять процессом расщепления урана, и были построены первые атомные электростанции. Наступила новая эра энергетики. Энергия мирового запаса урана сейчас превышает в десятки раз все существующие запасы энергии, рассмотренные ранее. Ситуация сейчас сложилась такая, что некоторые виды энергии начинают исчерпывать свой ресурс. В начале становления человечества и технологий люди потребляли энергию и не особо заботились о том, что будет далее. В настоящее время уже достаточно оскудели запасы угля и нефти. Возникли также всевозможные экологические проблемы. Один из возможных подходов решения этой задачи заключается в альтернативных источниках энергии, многие институты и компании занимаются этой проблемой. Одно из самых перспективных направлений в этой области является кумуляции энергии в слоистых системах. Данная проблема затрагивает различные области математики и физики. Поэтому для решения этой задачи важно развивать как теоретические подходы, так и экспериментальные. Однако эксперименты в данной области очень дороги, а порой и невозможны, так как с технической стороны еще не существует таких источников тяжелых ионов, которые смогли бы обеспечить необходимую степень кумуляции энергии в термоядерной мишени, хотя такие работы ведутся. Поэтому, на данном этапе математическое моделирование кумуляции энергии в слоистых системах оболочек является важнейшей задачей на пути к созданию термоядерных мишеней.
Работы по решению проблемы управляемого термоядерного синтеза ведутся в настоящее время в крупнейших лабораториях мира в таких стра-
нах, как: США, Франция, Япония, Германия и Россия. Наибольшие средства в настоящее время затрачиваются на стационарные установки с магнитным удержанием (ТОКАМАКи), лазерный и тяжелоионный термоядерный синтез с инерционным удержанием плазмы. Работы по лазерному и тяжелоионному термоядерному синтезу требуют создания установок с энергией выхода в несколько мегаджоулей и соответственно миллиардных затрат.
Тяжелоионный термоядерный синтез является наиболее перспективным направлением в данной области, так как энергия выхода в такой системе гораздо выше, чем в других. Для реализации концепции кумуляции энергии в слоистых системах необходимо разработать мишень и камеру, в которой будет проходить процесс обжатия мишени, а также установку на тяжелых ионах [1], способную вложить необходимую энергию в эту мишень.
Теоретической основой для реализации тяжелоионного термоядерного синтеза послужили работы по кумуляции энергии. Еще Релеєм и Гюгонио было показано, что используя класс автомодельных волн Римана, можно в процессе изоэнтропического сжатия плоского слоя политропного газа получить сколь угодно большую плотность газа [2]. Кумулятивные процессы и способы их реализации стали важнейшим направлением исследований в нашей стране с середины пятидесятых годов. Разнообразные типы кумуляции энергии были глубоко изучены академиком Е.И.Забабахиным, который привлек к этой проблеме большую группу теоретиков и экспериментаторов. Е.И.Забабахиным были рассмотрены различные случаи кумуляции энергии при фокусировке сферических и цилиндрических оболочек и ударных волн, дан анализ влияния вязкости, теплопроводности и сжимаемости на характер кумуляции. Были рассмотрены сходящиеся волны в веществах с фазовыми переходами и ударные электромагнитные волны, решены задачи о кумуляции энергии при несимметричной фокусировке и в не сходящейся ударной волне. Был открыт новый класс автомодельных кумулятивных течений в слоистых периодических системах с чередующимися слоями из легких и тяжелых ве-
ществ. В таких системах степень возрастания энергии на фронте ударной волны выше, чем в однородном веществе. Также в работах [3],[4] было отмечено, что процессы безударного сжатия газа являются энергетически более выгодными, так как не приводят к большому росту кинетической энергии и сильному разогреву вещества, что наблюдается при ударном сжатии. Продолжению исследований о безударном сжатии и кумуляции энергии в рамках модели идеального газа были посвящены работы А.Ф. Сидорова [5]- [8] и ряда других авторов, где рассматривались следующие задачи о безударном коллапсе:
движением непроницаемых подвижных сжимающих поршней сжать в точку массу газа, занимающую первоначально заданный объем, при этом получить неограниченно растущую плотность вещества и найти законы управления движением поршней;
получить в заданный момент времени требуемую степень сжатия исходной массы газа с наименьшими затратами энергии на движение сжимающих поршней, конфигурация которых известна в начальный момент времени.
Обзор данных работ по кумуляции энергии при безударном сжатии газов приведен в [9]. В [9] выделяется два периода исследований. На первом этапе изучаются автомодельные задачи, для цилиндрических и сферических течений, определяются степени кумуляции энергии для основных газодинамических величин. К примеру: в работах [10],[11] рассматриваются автомодельные режимы сжатия конечной массы плазмы. Также с использованием аналитических конструкций строятся соответствующие законы оптимального управления процессом с одной точкой переключения управления. Строятся не автомодельные режимы безударного сжатия плоских газовых слоев для получения локального роста плотности газа при конечных затратах энергии. Исследуются процессы безударного сжатия призм, тетраэдров. Идет поиск оптимальных форм мишеней, в которых при меньших энергозатратах до-
»
стигаются более высокие степени кумуляции плотности и энергии, чем при одномерном сферическом подходе. Рассматриваются вопросы устойчивости процесса сжатия при возмущении законов управления сжатием. Исследования первого этапа показали, что существуют формы мишеней, для которых при меньших энергозатратах можно достичь более высокие степени кумуляции плотности и энергии, чем при одномерном сферическом подходе. Все процессы первого этапа рассматривались в рамках модели идеального газа, без учета диссипативных процессов и других физических эффектов. Эти обстоятельства явились исходными для детального исследования многомерных процессов с целью конструирования мишеней новых нетрадиционных форм и учитывающих диссипативные процессы и другие существенные физические эффекты. Поэтому на следующей стадии исследований использовались численные методы. Первые расчеты сжатия призм, тетраэдров и конусообразных объемов газа проводились по методикам сквозного счета. Газодинамические расчеты по серийным методикам подтвердили полученные аналитические эффекты сверхкумуляции, а также показали, что неадаптированные специально для задач сверхсжатия методики позволяют решать достаточно надежно двумерные задачи до максимального сжатия в несколько тысяч раз и трехмерные в несколько сот раз. Был проведен анализ влияния лучистой теплопроводности для случая сжатия плоского слоя [12]. Была выполнена приближенная оценка степени кумуляции плотности с учетом теплопроводности, которая показала, что эта степень увеличивается в полтора раза по сравнению с аналитическим решением Релея - Гюгонио. Расчеты также показали, что среднее сжатие с учетом теплопроводности при одинаковых затратах вложенной энергии ухудшается, максимальное сжатие части слоя у поршня намного выше и кумуляция плотности увеличивается. Также отмечается, что в [13] были, приведены расчеты безударного сжатия газовой призмы и тела вращения, у которых есть неподвижная жесткая стенка. Расчеты проводились как без учета, так и с учетом горения. Результаты сравнивались с аналитическими
решениями, полученными в [7]. Было показано, что до начала энерговыделения энтропийная функция остается близкой к начальному значению, т.е. сжатие осуществляется в безударном режиме. Во всех рассмотренных выше случаях в сжимаемых объемах присутствовали неподвижные непроницаемые стенки, реализовать которые в экспериментах затруднительно из - за возник- новения зон больших давлений. Поэтому предпочтение было отдано изучению процессов неограниченного безударного сжатия замкнутых объемов газа специальных форм, на внешнюю границ которых действуют распределенное по времени и по поверхности давление. Исследования в этом направлении проводились как аналитически, так и численно. Были исследованы особенности нестационарных полей течений, возникающих при безударном сжатии осессиметричных конусообразных со сферическим сектором объемов полит-ропного газа. В работах [14],[13] под руководством А.В. Забродина впервые были рассчитаны безударные сжатия и горения двумерных замкнутых конструкций DT - газа. Горение газа при этом происходило с высоким энерговыделением. На основе этих построений численно конструировались мишени. Для реализации безударного сжатия DT - газа было найдено распределение энерговложения по направлениям. В [14],[13] также отмечалось, что в при заданном асимметричном энерговложении не разрушающуюся кумуляцию получить не удалось. В [9] также были выделены проблемы, которые по сей день остаются открытыми. Эти задачи заключаются в следующем: математическое проектирование физического эксперимента по реализации этого подхода к неодномерному безударному сжатию газов, поиск оптимальных путей реализации этого процесса при использовании различных физических полей для сжатия вещества. Отмечаются также определенные результаты, достигнутые при расчетном конструировании оболочечных слоистых микромишеней на основе реализации концепции безударного сжатия. В большинстве из работ кумуляция реализуется за счет геометрических особенностей картины течения. В [15] впервые был опубликован пример сконструирован-
ной кумуляции, когда ударная волна проходила по системе из чередующихся плоских слоев тяжелого и легкого веществ, толщины которых подбирались специальным образом.
Физика развития гидродинамических неустойчивостей, приводящих к нарушениям симметрии сжатия мишеней оказалась одной из наиболее важных, трудных и интересных проблем в инерциальном термоядерном синтезе. В ИПМ им. М.В. Келдыша [17] исследована неустойчивость Релей-Тейлора, возникающая при обжатии термоядерных мишеней. В работе [16] дан анализ экспериментальных и теоретических работ зарубежных и отечественных лабораторий в этой области. В частности, представлены результаты натурных и вычислительных экспериментов по изучению роста неустойчивостей в плоской и цилиндрической геометриях при действии лазерного или рентгеновского излучения на образцы. Эти данные связывают характеристики развивающихся возмущений с отклонениями от симметрии в начальных условиях, внесенными неоднородностями облучения или технологией изготовления мишеней. Стоит также отметить обзор теоретических работ по ИТС, проведенных в РФЯЦ-ВНИИТФ, где представлены результаты расчетов мишеней для ИТС, проведенных по одномерной программе ЭРА и комплексу двумерных программам ТИГР-ОМЕГА-ЗТ. Расчетно подтверждена возможность достижения коэффициентов усиления по энергии около 100 при быстром инициировании мишеней с помощью лазеров с энергией менее 500 кДж и длительностью импульса 20-50 псек.
В работах [18]-[21] были рассмотрены движения оболочечных слоистых систем при мгновенном энерговложении во внешние слои. Были построены последовательности асимптотик, описывающие основные закономерности движения систем при схождении к центру. Построенные последовательности приближенных решений для вычисления величины кумулирующейся энергии в слоистых системах не только описывают качественную зависимость от исходных параметров, но и с хорошим приближением дают количественные
значения. С использованием полученных результатов в [22], [23] была исследована возможность реализации безударного сжатия одномерных слоистых систем. Также было установлено, что для практического осуществления безударного сжатия газа можно построить слоистую оболочечную и реализовать энерговложение в один из внешних слоев так, чтобы фактически оболочечную систему можно было рассматривать как мишень, внутри которой будет происходить безударное сжатие [24]. Показано, что требуемые параметры по температуре и плотности при минимальном вложении энергии, необходимом для зажигания мишени, достижимы. Проведено сравнение с известной мишенью ИТИС [25],[26].
При обжатии термоядерных мишеней до высоких плотностей и энергий инициируется процесс их загорания. Возгорание мишени напрямую зависит от энергии, сконцентрированной в горючей смеси. В работах [27] - [30] построена математическая модель, реализующая уравнения одномерной радиационной гидродинамики с теплопроводностью и кинетикой термоядерного горения. Эта же модель используется в программном комплексе НЗТ для реализации кинетики термоядерного горения. При безударном сжатии температура горючей смеси недостаточно высока для возгорания, поэтому начиная с некоторого момента времени нарушают этот режим и энергию вкладывают по закону Q = Qmax- При этом возникающие ударные волны повышают температуру горючей смеси. В [31] приведены результаты расчетов с различными параметрами энерговложения.
На основе вышеперечисленных работ была разработана концепция реализации безударного сжатия [32], которая позволяет учитывать конструктивные характеристики установки, на которой планируется проводить эксперимент. Выделяются две стадии предложенной концепции: первая - теоретическая, целью которой является определение величины кумулирующейся энергии и выявления ее зависимости от параметров конкретной конструкции микромишени и способов энерговложения. Указанные зависимости позволяют предва-
рительно исследовать различные способы организации мишеней и отобрать наиболее подходящие из них по выходным параметрам кумуляции, что позволяет существенно сузить множество вариантов, которые следует рассматривать на последующей стадии. На следующей стадии численными методами исследуются мишени, отобранные на первой стадии, а именно те, которые в первом приближении удовлетворяют конструктивным характеристикам установки и выходным параметрам кумуляции энергии. Целью численных расчетов является подтверждение величины кумулирующейся энергии при обжатии мишени с использованием более точных математических моделей, а также всестороннее ее изучение. Такого рода задачи очень трудны и требуют для своего решения использования точных численных методов. В институте прикладной математики им. М.В. Келдыша под руководством А.В. Забродина разрабатывается многофункциональный программный комплекс НЗТ для расчета нестационарных двумерных уравнений газовой динамики с теплопроводностью на подвижных криволинейных сетках. Для возможности изучения влияния диссипативных процессов и других физических эффектов на кумуляцию энергии и прочих выходных параметров мишени необходимо было поставить соответствующую математическую модель. Математическая модель комплекса НЗТ изложена в [33], плазма представлена в виде трех компонент - ионов, электронов и фотонов, движущихся с общей скоростью, но имеющих разные плотность, давление и энергию. Данная модель позволяет учитывать теплопроводные и релаксационные процессы, проходящие в трехкомпонентной плазме в приближении серой материи.
Как уже отмечалось ранее, решение задач термоядерной физики требует применение как точных математических моделей, способных учесть основные физический процессы, происходящие в мишени, так и точных численных методов. Как правило при написании разностных схем для такого рода задач применяют метод расщепления по процессам [34],[35],[36]. Шаг по времени равномерно дробится на количество равное количеству выделенных процес-
сов. На первом шаге рассчитывается первый физический процесс, определяются неизвестные величины. Далее шаг за шагом рассчитываются остальные физические процессы, результаты которых рассматриваются как поправки к величинам найденным на предыдущем шаге. Использование данного метода позволяет реализовывать для разных физических процессов различные по своей идеологии разностные схемы. В программном комплексе НЗТ используется расщепление исходной математической модели на газодинамический и теплопроводный этапы.
Для расчета газодинамики существует огромное множество различных методов, каждый из которых обладает как положительными, так и отрицательными чертами. Можно выделить два основных подхода: первый - схемы, основанные на методе Годунова [45], [37] - [43]; второй - схемы сквозного счета [45], [44]. Рассмотрим каждый из подходов по отдельности. В схемах, основанных на методе Годунова, общая идея заключается в расчете задачи о распаде произвольного разрыва - задаче Римана, точно или приближенно. Результатом расчета являются промежуточные газодинамические величины. Они используются как величины промежуточного слоя по времени в аппроксимации исходной системы. Различия схем данного подхода заключаются в способе написания аппроксимации исходной системы, методов счета распадов разрыва, а также способе задания функции для газодинамических величин в пределах ячейки: постоянная, кусочно - линейная или кусочно - полиномиальная функция, TVD схемы реконструкции функций [46]. Для расчета газовой динамики в НЗТ используется обобщенный метод Годунова с выделением разрывов. Этот метод является первого порядка точности по пространству и времени, позволяет выделять основные разрывы. Использование разностных методов сквозного расчета газовой динамики без применения специальных мер не может в достаточной степени гарантировать сохранение энтропии и потенциальности движения. Более подробный обзор по схемам газовой динамики изложен в работах [44],[47].
При расчете сложных газодинамических задач возникает необходимость в построении на каждом временном слое адаптивной сетки. Простейшие алгоритмы построения сетки приведены в [45]. Более сложные методы описаны в следующей литературе: построение ортогональных разностных сеток посредством расчета конформных отображений [48]-[54]; итерационные методы построения сетки, с минимизацией некоторого функционала, отвечающего за геометрические особенности сетки [55]. В комплексе НЗТ используются итерационные методы построения сетки, а также методы описанные в [45]. При решении задачи о вытекании вещества с торцов мишени возникают весьма сложные геометрические конфигурации, что и приводит к необходимости использовать итерационные методы построения сетки. Функционалы, отвечающие за геометрические особенности сетки, имеют свободные параметры, управляя которыми возможно изменять внутреннюю структуру сетки. Данный произвол позволяет избегать большинства неприятных ситуаций с сеткой, таких как перехлестывание ячеек и т.п., однако иногда это и не помогает. Таким образом, проблема построения адаптивных сеток до сих пор является актуальной. Также следует отметить, что при решении задач в сложных с точки зрения геометрии областях возникает проблема построения регулярной сетки. Один из способов преодоления данной проблемы является метод разбиения исходной области на топологические четырехугольные подобла-сти(сіотаіп decomposition) [56]-[58]. При этом, в каждой из подобласти сетка строится регулярным образом вышеперечисленными методами. Данный метод позволяет сохранить регулярность сетки в расчетной области, но одновременно ставит проблему сшивки решения на границе двух подобластей.
Как уже отмечалось ранее, на теоретических этапах исследований термоядерных мишеней полагалось, что газ является идеальным, в реальности же плазма имеет весьма сложные уравнения состояния трехкомпонентной среды. Поэтому возник вопрос, как задавать эти уравнения состояния. Один из подходов заключает в сплайн - аппроксимации уравнений состояний [59]. В [33]
для трехкомпонентіюй среды: ионы, электроны, фотоны, в рамках модели комплекса НЗТ, выполнено локальная аппроксимация двучленными уравнениями состояния.
Для аппроксимации теплопроводных членов в задачах термоядерной динамики, как правило, используют консервативные схемы [60],[61]. На теплопроводном этапе комплекса НЗТ учитываются теплопроводные и релаксационные процессы. Для аппроксимации теплопроводных членов комплекса НЗТ в работах [62], [63] была построена консервативная разностная схемы на основе локальных среднеквадратичных приближений на криволинейных сетках. В работах [62],[63], [25] было проведено исследование динамики мишеней тяжелоионного термоядерного синтеза в приближении теплопроводной газодинамики. Также стоит отметить метод опорного оператора, позволяющего построить консервативную разностную схему для теплопроводных членов на криволинейной сетке [64],[65].
Таким образом, для каждого выделенного этапа были построены разностные схемы, каждая из них в отдельности обладает консервативностью, устойчивостью и аппроксимирует исходные интегральные уравнения [66]-[68]. Из теории разностных схем [35] известно, что порядок аппроксимации не расщепленной системы будет иметь суммарный порядок аппроксимации, а также обладать консервативностью, если каждая из схем обладала данным свойством. Известно, что явная схема Годунова [45] устойчива при выполнении условия Куранта. Система уравнений теплопроводного этапа НЗТ является нелинейной, поэтому для аппроксимации системы ее предварительно линеаризуют. В итерационном алгоритме H3T-LIM (локально итерационный метод) итерации ведутся по линеаризованной системе, количество итераций определяется по теплопроводному числу Куранту как \ Кт* (К ~ коэффициент теплопроводности). Последнее условие при К » 1 приводит к необходимости считать систему уравнений с большим числом итераций.
Для задач тяжелоионного термоядерного синтеза важен этап верифика-
ции решений, полученных численными методами. При этом существенно количество и качество различных вариантов. Первый способ - это расчет на последовательности сгущающихся сеток [6б]-[68]. Второй способ - это реализация нескольких методов аппроксимации исходной системы, что позволяет более широко судить о качестве получаемого решения, а также дает понять, какие эффекты вносятся численной схемой, а какие действительно имеют место. Третий способ - сравнение результатов с результатами полученными другими программными комплексами, основанных на общих математических моделях, например с программным комплексом СНД[69], разработанным под руководством Г.В. Долголевой для решения одномерных задач тяжелоионного термоядерного синтеза.
Теоретические модели играют важную роль на стадии предварительного проектирования, позволяют выявить диапазон допустимых параметров для дальнейшего их изучения на основе численных расчетов. Таким образом, разработка, развитие и верификация существующих программных комплексов, позволяющих рассчитывать широкий класс задач термоядерного синтеза, является наиболее актуальной задачей.
Диссертационная работа посвящена:
Совершенствованию методики НЗТ для расчета двумерных нестационарных течений теплопроводного газа в трехтемпературном приближении на криволинейных сетках;
Верификации методики и модели программного комплекса НЗТ;
Численному моделированию цилиндрических мишеней ИТИС с учетом термоядерного горения;
Настоящая работа продолжает цикл работ [18]-[33],[62],[63],[71],[72],[74], посвященных моделированию мишеней тяжелоионного термоядерного синтеза и является развитием программного комплекса НЗТ.
В главе 1 предлагается неявная разностная схема H3T-NLA (nonlinear algebra) для нелинейной системы уравнений теплопроводного этапа про-
1/
граммного комплекса НЗТ для счета на блочно-структурных сетках. Для решения СЛАУ линеаризованной системы НЗТ уравнений используется пакет программ Sparskit. Проведено сравнение результатов расчетов по модифицированной методике H3T-NLA и немодифицированной H3T-LIM (локально итерационный метод). Получено совпадение результатов. Преимущества схемы:
В H3T-LIM решается линеаризованная задача итерационным методом LIM, число итераций определяется НЗТ-числом Куранта; В H3T-NLA методике рассчитываются нелинейные уравнения, ведется контроль консервативности схемы; стратегия выбора шага по времени основана на контроле невязки;
На этапе загорания мишени H3T-LIM требует большого числа итераций. В H3T-NLA время счета шага практически постоянно; схема в целом экономичнее в 2 - 3 раза;
Реализованный алгоритм сшивки температур и потоков на границах ярусов обеспечивает консервативность при любых соотношениях длин смежных границ счетных блоков;
В Главе 2 рассматривается алгоритм построения фронтов ударных волн и контактных границ с искусственной вязкостью, алгоритм предложен С.К. Годуновым для расчета выделенных стационарных волн, проводятся тестовые расчеты. В продолжение исследований влияния открытых торцов мишени, на основе нового алгоритма движения контактных границ проведены расчеты по НЗТ методике для урановой мишени, которые показали, что есть существенно - одномерная зона. Новый алгоритм движения контактных границ позволил избежать вычислительных трудностей, свойственных старому алгоритму. Преимущества алгоритма:
Новый алгоритм, вносит контролируемую вязкость порядка 0(h) и га
сит нефизичные высокочастотные гармоники
В Главе 3 приведены результаты прикладных расчетов, целью которых является верификации методики и модели программного комплекса НЗТ:
На последовательности сгущающихся сеток в НЗТ рассчитана задача обжатия мишени с учетом термоядерных реакций в урановой оболочке (B.C. Имшенник, В.Т. Жуков), показано, что интегральные величины сходятся с 0(h). Проведен анализ влияния термоядерных реакций на параметры мишени;
На серии ID задач получено совпадение результатов с ID методикой SND ( Г.В. Долголева );
На задаче с энерговложением F—30, Q*=50 без учета горения получено совпадение результатов НЗТ и SND (объемный случай); В задаче с энерговложением F=30, Q*=50 (плоский случай) с учетом горения получено качественное совпадение результатов;
В 2D RZ - геометрии исследована формула безударного сжатия ( А.В. Забродин ), полученная для плоского случая. Показано, что она может быть применена для RZ - геометрии с энерговложением по объемным или плоским массам( последний вариант эффективнее)
Реализована программа оптимизации модели, она позволяет подобрать оптимальную толщину слоев для обеспечения более высокого энерговыделение в DT; показана практическая оптимальность выбранной конструкции мишени;
B.C. Имшенником предложена модель, учитывающая неравномерность излучения энергии фотонами по углам и частотам (2Т-1Р). В DT - слое вместо уравнения теплопроводности для фотонов решается уравнение переноса. На основе новой модели реализована методика Н2Т-1Р. Проведены расчеты для F=30, которые показали, что для ЗТ модели с объемным энерговложением точка загорания Q*=78.5, т.е. при меньших Q* горения нет. Расчеты Н2Т-1Р показали, что учет неравномерности излучения позволяет незначительно: понизить энерговложение в мишень
для ее загорания, при Q*=75 горение есть, что подтверждает трехтем-пературную модель, применяемую в комплексе НЗТ.
Основные результаты диссертационной работы сформулированы в заключении.
Результаты диссертационной работы опубликованы в работах [79],[80]-[82] и докладывались на пяти конференциях:
Н.М. Гиззаткулов. Алгоритм построения фронтов ударных волн и контактных границ с искусственной вязкостью. XV Всероссийская конференция Теоретические основы и конструирование численных алгоритмов для решения задач математической физики с приложением к многопроцессорным системам, посвященная памяти К.И. Бабенко (7-12 сентября,2004) Абрау-Дюрсо, Новороссийск;
А.В. Забродин, Н.М. Гиззаткулов. Разностная схема для двумерных нестационарных уравнений газовой динамики в трехтемпературном приближении. Тезисы докладов конференции молодых ученых "Ломоносов - 2004"МГУ. Москва. Апрель, 2004;
Н.М. Гиззаткулов, Численное моделирование двумерного нестационарного течения теплопроводного газа в трехтемпературном приближении, Тезисы докладов IX Всероссийского совещания по проблемам построения сеток и XIV Всероссийская конференция "Теоретические основы и конструирование численных алгоритмов для решения задач математической физики"(Дюрсо, 15 - 22 сентября 2002);
Н.М. Гиззаткулов, М.Ю. Заславский, А.Х. Пергамент. Регуляризован-ные алгоритмы решения уравнения теплопроводности на криволинейных сетках. Тезисы докладов конференции молодых ученых "Ломоносов - 2003"МГУ. Москва. Апрель, 2003;
Т.Т. Гарипов, Н.М. Гиззаткулов, А.Х. Пергамент, Эволюция напряженно - деформированного состояния пористых насыщенных сред в про-
цессе нефтедобычи, Повышение нефтеотдачи пластов. Освоение труд-ноизвлекаемых запасов нефти. Труды 12- го Европейского симпозиума "Повышение нефтеотдачи пластов"(Казань, 8-10 сентября 2003 года);
Диссертация и отдельные ее части докладывались на семинарах в Институте прикладной математики им. М.В. Келдыша РАН.
Автор выражает благодарность А.В. Забродину, В.Т. Жукову, С.К. Годунову, B.C. Имшеннику за предложенные задачи и конструктивные замечания, высказанные в процессе работы, В.Т. Жукову, Е.А. Забродиной и Л.А. Плинер за помощь при освоении комплекса НЗТ и при реализации новых вычислительных алгоритмов для него.
Система дифференциальных уравнений модели
Схема счета основана на декомпозиции (при необходимости) области решения на подобласти. Эта операция включает три этапа: геометрическое разбиение исходной расчетной области на подобласти (введением при необходимости дополнительных границ), движение границ, построение сетки в каждой подобласти. Разбиение исходной области на подобласти определяется, в частности, необходимостью ее представления набором "топологических четырехугольников". Эти "четырехугольники"должны в ходо расчета покрывать всю область искомого решения. Несколько подобластей могут объединяться в структуру, называемую ярусом, если при этом для объединения сохраняется топология четырехугольника, а сами подобласти составляют прямоугольную решетку. Ярус является основной вычислительной структурой: на каждом шаге по времени уравнения (2.1) - (2.5) решаются в ярусе с определенными краевыми условиями. На этапе счета теплопроводности несколько ярусов могут объединяться в счетные объекты называемые кластерами, в которых сквозным образом осуществляется расчет теплопереноса и релаксации.
Каждая из подобластей (счетная область) является либо физической областью, либо ее частью. Принадлежностями области являются ячейки, границы, физические параметры (параметры уравнения состояния, коэффици zoента теплопроводности и т.д.). В каждом ярусе мы строим четырехугольную криволинейную сетку. Особенностью метода является процедура учета взаимосвязи между ярусами, не предполагающая гладкости сетки через границы ярусов. В данном параграфе будет описана модификация теплового этапа расчета.
Система дифференциальных уравнений = diviKigiadT + CeiiTe + Qi, m = div(tfegradTe) + cei{Ti - Te) + Cef(T} - Ге4) + Qe, J dt дЕе{Те) дЕііТі) P P- = div(KfgrzdTf) + cef(T?f) + Qf, описывает теплопроводные и релаксационные процессы. Разностная аппроксимация имеет вид Ё{ = Д- + -Д(7;-) + -Сй(Ге-Т0 + "ф, Р Р Р Ёе = Ee + -A(fe) + -cei(fi-fe) + -cef(f}-f?) + -Qej (3.6) Р Р Р Р Ef = Ef + -A(f}) + -Cef(fi}) + -Qf. Иг г Здесь г - шаг по времени, А - разностная аппроксимация теплопроводных, дивергентных членов, конкретный вид которых не принципиален. Энергии на верхнем временном слое линеаризуются по температуре: ЕІ = ОІТІ + /ЗІ, Ёе = aefe + &, Ef = afff + /3/; (3.7) соответствующие формулы для коэффициентов ak и /3k (к = і, е, f) выписаны в [62]. Температуру Те4 можно линеаризовать двумя способами: либо f 4 = Zlmfj + Zln ГДЄ
Решая полученную систему, получим вектор температуры Т на момент времени t-\. Линейная система является плохим приближением исходной нелинейной системы, поэтому для того, чтобы гарантировать хорошую аппроксимацию, необходимо организовать итерационный процесс по методу Ньютона [68]. А именно, температуру (Те) линеаризуем согласно (3.9) в окрестности е виде Тк 1. Тогда система (3.6) с учетом (3.7) и линеаризации переписывается в oaff Выполняя к итераций, получим вектор температур Тк, где Т=Т.
На первом этапе казалось, что контроля количества итерации по невязке достаточно для решения нелинейной системы. Практика показала, что в задачах физики плазмы при прохождении фронтом тепловой волны границы раздела двух существенно разных веществ абсолютные невязки после Ньютоновских итераций порядка 10 3 — 10+1, при начальной абсолютной невязке порядка (108-1010). Возникает характерная картина - прогрев газа в режиме тепловой волны, фронт которой движется с конечной скоростью. Расчет такого режима затруднен тем, что температура существенно негладкая около точки фронта. Поэтому необходимо было ввести алгоритм, который вел бы контроль шага по времени от величины невязки
Алгоритм переинтерполяции граничной скорости из полуцелых точек в целые
Выполняя к итераций, получим вектор температур Тк, где Т=Т.
На первом этапе казалось, что контроля количества итерации по невязке достаточно для решения нелинейной системы. Практика показала, что в задачах физики плазмы при прохождении фронтом тепловой волны границы раздела двух существенно разных веществ абсолютные невязки после Ньютоновских итераций порядка 10 3 — 10+1, при начальной абсолютной невязке порядка (108-1010). Возникает характерная картина - прогрев газа в режиме тепловой волны, фронт которой движется с конечной скоростью. Расчет такого режима затруднен тем, что температура существенно негладкая около точки фронта. Поэтому необходимо было ввести алгоритм, который вел бы контроль шага по времени от величины невязки. Известно, что метод Ньютона решения уравнения f{x) = 0 сходится квадратично ([93]) при условиях: rk C-\CrQf-\\fxx{x)\\ С2; \\f-\x)\\ Си С = C2Ch
Теперь можно более точно указать, насколько хорошим должно быть начальное приближение яо? чтобы процесс сходился. Очевидно, для этого достаточно выполнения неравенства: . С/(х) q 1.
В нашем случае для нелинейного уравнения не определялись константы С\ и Сг, а предполагалось, что линеаризация (3.8) позволяет выбрать начальное приближение таким образом, чтобы процесс решения нелинейного уравнения сходился. В случае плохого начального приближения инициировался процесс мельчения шага по времени для того, чтобы попасть в окрестность решения нелинейного уравнения. Исходя из вышесказанного, был предложен следующий алгоритм. На каждой Ньютоновской итерации рассчитываются невязки Г{, re, rf нелинейной системы (3.6), в норме Z/2 по счетному пространству. Итерации по нелинейности выполняются до тех пор, пока одна из невязок убывает в заданное число раз (один, два порядка или более) или не достигается величина птах = 10(максимальное количество Ньютоновских итераций). В случае, когда одна из невязок после итерации отличается от средней величины невязки, выработанной за предыдущие сто шагов, более чем в заданное число раз(103 — 104), то шаг по времени г уменьшается вдвое, инициируется процедура возврата на перерасчет газодинамического этапа с новым шагом по времени, процесс повторяется до тех пор, пока невязка не достигает средней величины. В [93] приведен ряд алгоритмов, позволяющих оптимальнее вырабатывать шаг приращения временной переменной. На данном этапе реализован самый простой из всех вариантов для проведения методических расчетов и верификации численной схемы.
Для решения системы линейных алгебраических уравнений с разреженной матрицей использовался программный комплекс Sparskit [75], использующий метод обобщенной минимальной невязки(СМ11Е8). GMRES считается одним из самых эффективных методов решения несимметричных линейных систем. В пакете использовались следующие настройки: є = 2.22 10 20 (абсолютный толеранс), rnaxits — 200, гт — 15, где im - размерность подпространства Крылова, maxits - максимальное количество итераций в методе обобщенных невязок. В процессе расчетов Sparskit зарекомендовал себя как надежный алгебраический решатель. К тому же, данный пакет имеет параллельный вариант Sparslib. О пакете Sparslib и его возможностях написано в [75].
На примере двух ярусов рассмотрим алгоритм сшивки решения теплопроводного этапа программного комплекса НЗТ. На рис. 1.1 изображена схема расчета теплового потока на ребре границы яруса 1; кружочками и точками обозначаются тепловые потоки на граничных ребрах ярусов 1 и 2 соответственно, направленные по нормали к ребрам. Приведем формулы для расчета теплового потока в точке 5 рис. 1.1. В точке 5 проводится нормаль, на которую проецируются температурное поле ярусов 1 и 2, с использованием линейной интерполяции:
Апробация алгоритма построения фронтов ударных волн и контактных границ с искусственной вязкостью на задаче о вытекании вещества через торцы мишени
По поставленной задаче были проведены расчеты, границы 4,7,10 передвигались по новому алгоритму параграфа 1. Вязкость подбиралась таким образом, чтобы влияние сглаживания сказывалось лишь на высокочастотных гармониках разложения в ряд Фурье нормальной компоненты скорости вдоль оси Y, в данном расчете параметр S = 0.1.
На рис.2.32 - рис. 2.35 приведены профили газодинамических величин на момент времени t — 102.2, которые показывают образование разгрузки с торцов мишени. На рисунках рис.2.36 - рис. 2.39 приведены профили газодинамических величин на момент максимального обжатия мишени t = 102.74, они демонстрируют сохранение существенной зоны, которую не затрагивают разгрузки с торцов мишени. На рис. 2.40, 2.41 для t = 102.77 приведена фаза разлета мишени. На рис. 2.40 видно, что скорость одномерной зоны мишени положительна и равна 87, в то время, как в области разгрузки скорости направлены в центр и равны -1, поэтому в окрестности точки X = 35 возникает повышенное давление рис. 2.41. На рис. 2.40 приведена сетка, проходящая через центры ячеек расчетной сетки; сложности счета по предыдущему алгоритму возникали на момент разворота границы 10 в окрестности точки X = 35. Использование нового алгоритма позволило избежать таких трудностей.
В ходе численных экспериментов с новым алгоритмом было выявлено, что параметр 8 следует выбирать основываясь на разложении нормальной компоненты скорости в ряд Фурье вдоль фронта волны таким образом, чтобы вносимая вязкость гасила лишь нефизичные высокочастотные гармоники. Завышение коэффициента 8 может привести к изменению решения в области.
Рассмотрен алгоритм построения фронтов ударных волн и контактных границ с искусственной вязкостью, алгоритм предложен С.К. Годуновым для расчета выделенных стационарных волн, проведены тестовые расчеты по обтеканию сферы. На основе нового алгоритма движения контактных границ проведены расчеты по НЗТ методике для урановой мишени, которые показали, что есть существенно - одномерная зона. Новый алгоритм позволил избежать вычислительных трудностей, свойственных старому алгоритму.
В данной главе приводятся результаты расчетов некоторых задач, на которых были опробованы методы, разработанные в предыдущих главах. В первом параграфе ставится задача об обжатии термоядерной мишени, где в качестве "пушера"(слоя, сжимающего дейтерий - тритиевую смесь) выступает уран -238 [95]. Расчеты ведутся по алгоритму, предложенному в главе 1. Предполагается, что "делящаяся" урановая оболочка будет способствовать увеличению энерговыделения рабочей DT - области, сдерживая ее разлет. Задачи по повышению энерговыделения в DT - области являются наиболее актуальными на сегодняшний день, так как это позволит уменьшить энергию, которую следует вложить в мишень для достижения достаточного КПД (здесь КПД есть величина, равная отношению выделенной энергии от горения топлива к вложенной энергии ), также это позволит уменьшить требования к самой энергетической установке, основным параметром которой является величина Qmax пиковое обострение энерговложения, которое может выдать установка. Задача рассчитывалась на последовательности сгущающихся сеток. Демонстрируется, что интегральные величины сходятся с 0(h). Проводится анализ влияния термоядерных реакций на параметры мишени.
В параграфе 2 приводится сравнение двух методик расчета: СНД - движения двухтемпературного излучающего газа в одномерном случае ([69]), (ВНИИЭФ, рук. Г.В. Долголева), и НЗТ - движение трехтемпературного газа в двумерном случае ([62],[63],[72]), ( ИПМ, рук. А.В. Забродин), на одной задаче. Методики СНД и НЗТ основаны на единой,общепризнанной физиками, физической модели, а именно, газовая динамика с трехтемпературной теплопроводностью. Отличие методик заключается в реализации конечно -разностного счета газовой динамики с трехтемпературной теплопроводностью, а также в уравнениях кинетических термоядерных реакций.
Сравнение двумерной методики НЗТ с одномерной методикой СНД
Сравнение проводилось с момента времени t = 50 до момента разлета мишени. Момент t = 50 - начало движения урановой оболочки ( см. рис. 3.35), был выбран для того, чтобы отследить профили сравниваемых величин до момента t = 102.2 выхода возмущений на ось симметрии. На рисунках рис. 3.5 - рис. 3.33 приведены зависимости плотности, электронной темпера туры и скорости от радиуса в DT - слое термоядерной мишени на различные моменты времени (штрихпунктирной линией обозначается результат расчета по методике СНД, непрерывной линией - результат расчета по методике НЗТ на близкие моменты времени). В методике СНД при расчетах в DT слое использовалось 50 лагранжевых точек, в НЗТ - 160 точек, распределенных по геометрической прогрессии от центра мишени с показателем Л = 0.98. До момента t = 102.2, выхода возмущений на центр (рис. 3.26 - рис. 3.28) можно наблюдать практическое совпадение результатов расчетов по обеим методикам. Показанные на рис. 3.5 - рис. 3.28 профили температур, плотностей и скоростей практически совпадают. После выхода возмущений на центр (рис. 3.29 - рис. 3.31) есть небольшой сдвиг по времени, равный 0.03 единицы, это связано с тем, что газовая динамика в методике СНД считается со вторым порядком точности, а в методике НЗТ с первым, что в конечном счете влияет на скорость движения границ. На рисунках 3.32, 3.33 приведены профили температур и плотностей на момент максимального обжатия (102.6 для НЗТ и 102.63 для СНД). На рис. 3.35 приведена сравнительная RT диаграмма для обеих методик, которая демонстрирует совпадение движения всех слоев мишени.
Замечание. В процессе сравнения оказалось, что в методике СНД в повседневных расчетах использовалась отличная от методики НЗТ формула для энерговложения (более эффективный вариант), а именно были различны константы то,ші,Ш2,тз(см. Глава 1), задающие вид профиля энерговложения. В методике СНД: гао = Ро(го); mi = /oi(ri-ro); (2.1) т2 = Р2ІГ2 - n); rn3 = рз(г3 - г2). В методике НЗТ: (2.2) гао = Ро(го); mi = pi(rf-rg); т2 = Р2(г2-г1);т3 = рз{г%-г$).
Здесь ГІ, і = 0..3 - внешний радиус г +l слоя мишени, а рі, і = 0..3 - плотность вещества г + 1 слоя. Поэтому при сравнении использовался вариант вложения НЗТ. Далее будем называть эти энерговложения плоским и объемным, соответственно.
А.В. Забродиным аналитическим путем была получена формула энерговложения в случае плоской геометрии, для безударного сжатия мишени [9]. Реальные расчеты проводятся в 2D или ID RZ - геометрии, где аналитический вывод формулы энерговложения труден, поэтому применяется или формула, полученная для плоского случая или она слегка модифицируется и используется объемное энерговложение.
Для сравнения двух типов энерговложений были проведены одномерные расчеты в ID RZ геометрии для поставленной в параграфе 1.1 задачи, отличие заключалось лишь в том, что в расчетах использовались следующие параметры для энерговложения: Q = 50, F = 30, учет горения был включен. Вычисления проводились по двум методикам СНД и НЗТ с плоским и объемным энерговложениями. Результаты расчетов приведены в таблице 3. Проведенные вычисления: подтвердили качественное совпадение результатов расчетов СНД и НЗТ с учетом горения, т.е. мишень горела при плоском энерговложении и не горела при объемном; показали, что формулы для плоского и объемного энерговложения могут быть применены для 2D и ID RZ - геометрии. Отличие коэффициентов усиления в таблице 3 объясняется разными моделями кинетики термоядерных реакций.
Изучение газодинамических профилей для объемного и плоского энерговложений позволило ответить на вопрос: сохраняется ли безударное сжатие при применении, формул предложенных А.В. Забродиным для плоского случая, в ID и 2D RZ - геометрии. На рис. 3.36 - 3.39 представлены распределение скорости, плотности по радиусу DT - слоя мишени на различные моменты времени из интервала t = 0..99.97 (t = 99.97 выход из безударного режима сжатия мишени) для плоского и объемного энерговложений. Рис. 3.36 - 3.39 демонстрируют наличие слабых ударных волн в DT - слое, что говорит о том, что режим безударного сжатия нарушается. Проведенные расчеты (см. Таблица 3), показали, что формулы для объемного и плоского энерговложения применимы в случае ID и 2D RZ - геометрии. Формула плоского энерговложения более эффективна. Ее эффективность объясняется тем, что в DT слое мишени больший промежуток времени сохраняется высокая плотность и температура, так как урановая оболочка мишени движется медленнее (см. рис. 3.40). В случае объемного энерговложения достигаются более высокие плотности и температуры в DT - слое мишени, но общая длительность удержания топлива, при высокой температуре и плотности, меньше, поэтому DT - слой не успевает разгореться. На рис. 3.42, 3.41 приведены профили скоростей на момент максимального обжатия мишени. Рис. 3.41 демонстрирует торможение урановой оболочки на момент максимального обжатия для объемного энерговложения, в плоском случае урановая оболочка движется согласовано с топливом (рис. 3.42).
