Вычислительные технологии
Том 17, № 2, 2012
Один подход к численному моделированию безударного сильного сжатия специальных объемов газа
А. В. Рощупкин
Уральский государственный университет путей сообщения, Екатеринбург, Россия
e-mail: [email protected]
Рассмотрено численное моделирование процессов безударного сжатия идеального политропного газа до любой заранее заданной плотности. В расчетах, проведенных как в обратном, так и в прямом направлении изменения времени, получена степень сжатия газа в 103-104 раз.
Ключевые слова: газовая динамика, сильное безударное сжатие.
Введение
Решение задачи о безударном сильном сжатии газа рассматривается в связи с реализацией процессов управляемого лазерного термоядерного синтеза. Для получения реакций синтеза ядер вещества этому веществу требуется придать необходимые плотность и температуру. В настоящее время сформировался подход, согласно которому сжимаемое вещество представляет собой малую мишень, а энергия для сжатия и разогрева вкладывается лазерными установками.
Экспериментальные исследования в этом направлении масштабно ведутся, например, в Ливерморской национальной лаборатории им. Лоуренса, Ливермор, США [1]. В работе [2] предложена схема мишеней оболочечной структуры для лазерного управляемого синтеза. На рис. 1 схематично показа конфигурация мишеней. Стрелками обозначены направления вложения энергии.
В левой части рис. 1 изображена мишень цилиндрической формы, имеющая внешний прочный кожух и состоящая из среды, в которую будет вкладываться энергия, и среды, сжимаемой под воздействием расширения первой среды, разделенных выпуклой в сторону сжимаемой среды "треугольной" оболочкой. В качестве сжимающей среды может быть использован бериллий, в качестве сжимаемой — дейтерий-тритий, в качестве разделяющей оболочки — золото. Энергия в мишень вкладывается с торца. В правой части рисунка — мишень сферической формы, имеющая те же составные части. Сжимаемая и сжимающая среды разделены выпуклой в сторону сжимаемой среды оболочкой, имеющей форму "тетраэдра". Энергия в мишень вкладывается через четыре отверстия, расположенных против центров граней тетраэдра.
Отметим, что мишени другой оболочечной структуры рассмотрены в работе [3].
Предполагается, что если в предложенных в [2] мишенях разделитель сжимаемой и сжимающей сред будет двигаться в соответствии с решением задачи о безударном
* Работа выполнена при финансовой поддержке РФФИ (грант № 08-01-00052).
2В | 31)
Рис. 1. Мишень для лазерного термоядерного синтеза [2]
сильном сжатии, то сжимаемую среду удастся перевести в состояние с плотностью, необходимой для начала термоядерной реакции, без возникновения ударных волн.
В настоящей работе представлены два решения системы уравнений газовой динамики, описывающие процессы безударного сильного сжатия идеального газа до любой заранее заданной конечной плотности. Под термином "безударное" понимается, что все описываемые здесь течения газа отделяются друг от друга только слабыми разрывами, но не ударными волнами.
В пространстве физических переменных г, £ будем рассматривать одномерные из-энтропические течения, возникающие в политропном газе с уравнением состояния р = р7 /7, которые являются решениями системы уравнений газовой динамики
(т — 1) ( vu \
С + исг +----си +--=0,
2 2 У г) (1)
и + 7-т ссг + ииг = 0,
(7 — 1)
\у+1
где £ — время; г = 4/ ^ X > 0 — пространственная переменная; с = р(т-1)/2 — скорость
звука; 7 — показатель адиабаты; р — давление; V — параметр геометрии: V = 0 в случае плоских, V =1 в случае цилиндрических, V = 2 в случае сферических течений.
Задача 1. Пусть в начальный момент времени £ = £0 газ является покоящейся (и = 0) однородной средой с плотностью р =1. Это состояние назовем состоянием 1. Состояние газа в момент времени £ = > £0 с плотностью р, равной некоторому р* > 1, назовем состоянием 2. Требуется найти течения газа, возникающие при безударном переходе газового слоя из состояния 1 в состояние 2 под воздействием одного поршня.
Решение данной задачи в случае плоских течений впервые было предложено Мизе-сом, поэтому далее будем называть ее "конфигурацией Мизеса". В случаях одномерных сферически и цилиндрически симметричных течений, а также пространственных течений эта задача решена в [4]. Численное решение в сторону, обратную направлению изменения времени, построено в [5].
Задача 2. Пусть в начальный момент времени £ = £0 газ является покоящейся (и = 0) однородной средой с плотностью р =1. Это состояние назовем состоянием 1. В момент времени £ = > £0 газ однороден, его плотность р равна некоторому р* > 1; при этом газ покоится. Данное состояние газа назовем состоянием 2. Требуется найти течения газа, возникающие при безударном переходе одномерного газового слоя из состояния 1 в состояние 2 под воздействием двух поршней, один из которых покоится.
Задача 2 предложена в ряде работ, например, [6-10]. Далее будем называть ее задачей о сжатии "из покоя в покой".
В работах [11-13] доказано существование решения задачи о сжатии "из покоя в покой" в случае, когда газ сжимается двумя поршнями, один из которых покоится в точке с г = г0 > 0, а второй сжимает газ снаружи, т. е. координата второго поршня все время больше г0. Там же указано, что поскольку координата неподвижного поршня го строго положительна, то факт существования решения задачи о сжатии "из покоя в покой" имеет место и при сжатии газового слоя изнутри.
Численное решение задачи о сжатии "из покоя в покой" в сторону, обратную направлению изменения времени, построено в [14].
Итак, имеются решения двух задач о безударном сильном сжатии, построенных при изменении времени от к ¿0 < используя которые, можно получить траектории подвижных сжимаемых поршней, реализующих сжатие до требуемых величин плотности газа. Применение этих траекторий непосредственно к мишеням, предложенным в работе [2], — задача трудоемкая, поскольку требует двумерных и трехмерных нестационарных расчетов. Настоящая статья является началом ее реализации.
Рассмотрим численные решения системы уравнений газовой динамики, описывающие течения, возникающие при движении сжимающего поршня, траектория которого построена по алгоритмам из [5, 14], но в прямом направлении изменения времени.
1. Постановка задачи
Из рис. 1 видно, что сжимающие стенки выгнуты в сторону сжимаемой среды, т. е. речь идет о сжатии газа изнутри, поэтому далее будет рассматриваться только сжатие изнутри.
Алгоритмы, описанные в [5, 14] и реализованные с учетом нужного направления сжатия (снаружи или изнутри), в качестве одного из результатов дают траекторию контактной характеристики, через которую газ не течет. Данная поверхность трактуется как траектория сжимающего поршня, реализующего состояние 2 из задач 1 и 2. Траектория поршня строится в направлении, обратном изменению времени. По окончании построения ее координаты изменяются так, чтобы финальная точка переходила в точку с координатами г0 = 1, ¿0 = 0.
Теперь можно поставить задачу о поршне: пусть в начальный момент времени ¿0 = 0 в покоящийся однородный газ с плотностью р0 = 1 в точке с координатой г0 = 1 начинают плавно вводить поршень, заданный координатами точек с известными в них скоростями. В точке ¿0, гп, где координата гп определяется массой сжимаемого газа, в задачах 1 и 2 ставится неподвижная стенка, т. е. скорость газа равна нулю. Требуется численно построить течения газа, возникающие при движении заданного поршня.
2. Алгоритм расчета
Алгоритмом расчета является конечно-разностный метод, основанный на известном методе характеристик с пересчетом (см., например, [15, 16]). Решение строится в виде характеристической сетки в сторону прямого изменения времени.
При начале перемещения поршня вправо начинает двигаться звуковая характеристика С+, отделяющая покоящийся газ от уже сжатого (рис. 2). В задаче 1 данная
характеристика достигает неподвижной стенки в момент сильного сжатия £ = , а в задаче 2 она отражается от неподвижной стенки и в момент сильного сжатия должна вернуться на поршень. В силу этого различия в поведении характеристики С+ численные решения рассматриваемых задач различаются.
Рассмотрим алгоритм решения задачи 1 (см. рис. 2, а). Поскольку параметры несжатого газа известны полностью, то счет начинаем от характеристики С+, отделяющей покоящийся газ от уже сжатого, в сторону траектории поршня. Для этого С+ делится с некоторым шагом по времени Д£ точками с1, с2,..., с^,..., сп. Расчетная сетка строится слоями, где под слоем понимается численно построенная характеристика С-, выходящая из некоторой точки С+.
Из первой точки С1 до пересечения с траекторией поршня выпускается характеристика С-. Точка пересечения определяется пересечением характеристики С- и отрезка, соединяющего две соседние точки траектории поршня. Расчет инвариант Римана, а значит и газодинамических параметров, в ней ведется по формулам
Ьр = ¿1 + V (¿р - ¿1 )(7 - 1)
К - ¿1
8г1 :
К
2пр -¿р,
где Кр,Ьр — инварианты Римана в точке пересечения с траекторией поршня, ир — скорость газа на поршне, полученная интерполяцией по крайним точкам отрезка его траектории.
Из второй точки с2 выпускается характеристика С- до пересечения с характеристикой С + , выходящей из точки пересечения предыдущего слоя и траектории поршня. Расчет газодинамических параметров в точке пересечения ведется по стандартным формулам метода характеристик с пересчетом. Из этой новой точки до пересечения с поршнем выпускается характеристика С-. Газодинамические параметры в точке пересечения с поршнем рассчитываются по приведенным выше формулам.
Построение дальнейших слоев ведется до тех пор, пока характеристики С-, выходящие с С+ , пересекаются с траекторией поршня.
Рассмотрим алгоритм решения задачи 2 (см. рис. 2, б). В этом случае характеристическая сетка строится в два этапа. Первый этап: построение сетки и расчет газодинамических параметров в ее узлах в области течения, примыкающего к покоящемуся несжатому газу, — до точки С. Сетка на данном этапе строится аналогично алгоритму решения задачи 1. Второй этап: расчет сетки от точки С и выше вдоль стенки С В в направлении увеличения времени.
А А В
/с;
/АС"/ КЛ2
—►
г г г
' п ' 'о
г„ г
а
Рис. 2. Построение сетки для задач 1 (а) и 2 (б)
Пусть С — точка отражения характеристики С+ от неподвижного поршня в точке гп. Характеристика Свыходящая из С, строится по алгоритму первого этапа. Пусть О — первая точка, построенная на этой характеристике. Тогда для начала построения следующего слоя — характеристики С-, выходящей со стенки СВ, найдем точку пересечения характеристики С + , выходящей из точки О, и стенки СВ.
Зная уравнение характеристики С + , выходящей из точки О, находится точка ее пересечения со стенкой СВ, а инварианты Римана с учетом условия непротекания на поршне рассчитываются по формулам
—2 — Т2
—ь = —п - V(Ьь - 1В)(7 - 1)—^-£, Т = -ЯЬ;
где —ь, Ть — инварианты Римана в точке пересечения со стенкой.
При программной реализации рассмотренных алгоритмов расчета необходимо следить за размерами получаемых ячеек сетки. Если — некоторая заданная величина, то допустимыми ячейками сетки являются такие, размер сторон которых составляет не более 1.5Д£ и не менее 0.5Д£. В случае нарушения первого условия вводятся дополнительные точки, второго — точки удаляются.
Таким образом, имеются два параметра, влияющие на точность расчета: ДЬ — шаг разбиения характеристики С+, и Д£ — размер сетки.
3. Результаты расчетов
В качестве начальных данных принимался закон движения поршня, вычисленный по методике [5, 14], с относительной погрешностью масс сжатого и несжатого газа менее 1%.
При решении задачи 1 (конфигурация Мизеса) никаких особенностей выявлено не было. Как и в расчете при обратном направлении изменения времени, в данном случае плотность газа на поршне растет до ожидаемых величин, составляющих 103-104 начальной плотности. Характеристики одного семейства при построении сетки не пересекаются.
Иная картина наблюдается при численном решении задачи 2 (сжатие "из покоя в покой"). В качестве примера рассмотрим случай сжатия изнутри при 7 = 5/3, V = 1, т =1, где т — масса сжимаемого газа, для которого в табл. 1 приведены варианты построения траектории сжимающего поршня по алгоритму [14]. Здесь р* — конечная плотность сжатия, Д£, ДЬ, п — параметры, отвечающие за точность расчета (первые
Таблица 1. Варианты построения траектории поршня
Номер р* АС Аt п 5ш
варианта
1 1000 0.00001 0.00001 10 000 0.6408
2 1000 0.000005 0.000005 50 000 0.1766
3 1000 0.000005 0.000005 100 000 0.1427
4 1000 0.000001 0.000001 100 000 0.1068
5 1000 0.000001 0.000001 500 000 0.0159
6 1000 0.000001 0.000001 1000 000 0.0109
Р'-
58
56
54
52
Рис. 3. Временной срез значения плотности: а — вариант 1, б — вариант 6; кривая 1 для обоих вариантов взята из расчета в обратном направлении изменения времени
два имеют тот же смысл, что и в разделе 2, а п — число точек разбиение отрезка [1, с*]), 8 т — относительная погрешность масс сжатого и несжатого газа.
На втором этапе численного построения решения задачи 2 для всех вариантов (см. табл. 1) наблюдались пересечения характеристик одного семейства, т. е. возникали ударные волны. Чем больше была погрешность построения траектории сжимающего поршня, тем раньше появлялась ударная волна. В варианте 1 ударная волна возникала на первом расчетном слое второго этапа (рис. 3, а, кривая 2), а в варианте 6 в тот же момент времени ударной волны не наблюдалось (рис. 3, б, кривая 2) — она смещалась выше по времени.
Полученные волны трактуются как ударные волны слабой интенсивности, поэтому течения и далее считаются изэнтропическими [17], что оправдывает применение метода характеристик без дополнительных модификаций после их возникновения. Отметим, что после ударной волны в расчете на этом слое не возникает так называемого эффекта типа "пила" (последовательные многочисленные осцилляции).
В табл. 2 для шести вариантов (см. табл. 1) приведены максимальные значения плотности газа на стенке С В, до которой удалось досчитать, и также величина параметра, характеризующего время возникновения ударной волны. Здесь р** — максимальная плотность, полученная на стенке С В, — ¿5 — разность финального времени и времени возникновения ударной волны.
Таблица 2. Максимальные значения плотности
Номер АС А* * * ** ~*в
варианта
1 0.0001 0.00001 999.60 0.0741347 0.0030718
2 0.0001 0.00001 999.70 0.0740293 0.0033635
3 0.0001 0.00001 999.38 0.0739519 0.0027098
4 0.0001 0.00001 998.30 0.0739132 0.0027493
5 0.0001 0.00001 999.30 0.0739003 0.0000141
6 0.0001 0.00001 999.30 0.0738961 0.0000098
Представляется, что при расчете задачи 2 имеются следующие причины возникновения ударной волны. Во-первых, если при численном построении решения задачи 1 аналог центрированной волны Римана сходится в точке B (см. рис. 2), то в случае задачи 2 — в точке A. Поэтому при построении решения задачи 1 нет необходимости в приближении к точке B на критические расстояния, тогда как в задаче 2 для того, чтобы максимально возможно продвинуться по стенке CB вверх, необходимо рассчитывать все семейство центрированных характеристик C-, выходящих с этой стенки. Во-вторых, возникновение ударных волн в различные промежутки времени в зависимости от точности задания траектории поршня доказывает большую неустойчивость решения задачи 2 по сравнению с решением задачи 1. Ударные волны (пересечения характеристик одного семейства) для задачи 1 возникали при задании траектории поршня, построенного с погрешностью более 10 %.
Таким образом, на основании проведенных исследований можно сделать следующие выводы.
1. В одномерном случае предложена методика, позволяющая в прямом направлении изменения времени построить течения, возникающие в газе при движении сжимающего поршня, реализующего переход из состояния 1 в состояние 2 при сжатии газа изнутри.
2. Выявлены особенности численного построения течений. Предположено, что возникновение этих особенностей связано с недостаточно точным воспроизведением траектории сжимающего поршня.
3. Используя закон движения сжимающего поршня, в численных расчетах получена плотность сжатого газа, составляющая 103-104 его начальной плотности.
Автор выражает благодарность проф. С.П. Баутину за плодотворное обсуждение результатов работы.
Список литературы
[1] Süter L.J., Lindl J., Atherton J. и др. Успехи НИФ в области зажигания // Тезисы Междунар. конф. "X Забабахинские научные чтения". Снежинск: РФЯЦ ВНИИТФ, 2010. С. 5.
[2] Баутин С.П. Об одной конструкции мишеней для управляемого термоядерного синтеза // Там же. С. 30-31.
[3] ДолголЕВА Г.В., Забродин А.В. Кумуляция энергии в слоистых системах и реализация безударного сжатия. М.: Физматлит, 2004.
[4] Баутин С.П. Математическая теория безударного сильного сжатия идеального газа. Новосибирск: Наука, 2007.
[5] Николаев Ю.В. О численном решении задачи безударного сильного сжатия одномерных слоев газа // Вычисл. технологии. 2001. Т. 6, № 2. C. 104-109.
[6] Крайко А.Н. О свободном нестационарном расширении идеального газа // Изв. РАН. МЖГ. 1993. Т. 4. C. 155-163.
[7] Крайко А.Н. Вариационная задача об одномерном изэнтропическом сжатии идеального газа // Прикл. математика и механика. 1993. Т. 57, вып. 5. C. 35-51.
[8] Крайко А.Н. Асимптотические закономерности нестационарного расширения идеального газа в пустоту // Там же. 1994. Т. 58, вып. 4. C. 70-80.
[9] Крайко А.Н. О неограниченной кумуляции при одномерном нестационарном сжатии идеального газа // Там же. 1996. Т. 60, вып. 6. C. 1000-1007.
[10] Крайко А.Н., Тилляева Н.И. Автомодельное сжатие идеального газа плоским, цилиндрическим или сферическим поршнем // Теплофизика высоких температур. 1998. Т. 36, № 1. C. 120-128.
[11] Баутин С.П. О существовании решений задачи А.Н. Крайко // Прикл. механика и техн. физика. 2000. Т. 41, № 3. С. 48-55.
[12] Баутин С.П. О возможности изэнтропического перехода от однородного покоя в другое однородное покоящееся состояние идеального газа // Докл. АН. 1998. Т. 362, № 5. С. 621-624.
[13] Баутин С.П. Математическое исследование безударного сжатия газа // Успехи механики. 2002. Т. 1, № 2. С. 3-36.
[14] Николаев Ю.В. Численное решение задачи А.Н. Крайко // Вычисл. технологии. 2005. Т. 10, № 1. С. 90-102.
[15] Жуков А.И. Применение метода характеристик к численному решению одномерных задач газовой динамики // Труды математического института. 1960. Т. 58.
[16] Рождественский Б.Л., Яненко Н.Н. Системы квазилинейных уравнений и их приложения к газовой динамике. М.: Наука, 1968.
[17] Зельдович Я.Б., Райзер Ю.П. Физика ударных волн и высотемпературных гидродинамических явлений. М.: Наука, 1966.
Поступила в 'редакцию 28 июля 2011 г., с доработки — 26 сентября 2011 г.