Комбинированный метод расчёта безударного сильного сжатия одномерных слоёв газа в конфигурации Р. Мизеса Текст научной статьи по специальности «Математика»

структуры и моделирование 2017. №1(41). С. 93-102

УДК 519.6:533

КОМБИНИРОВАННЫЙ МЕТОД РАСЧЁТА БЕЗУДАРНОГО СИЛЬНОГО СЖАТИЯ ОДНОМЕРНЫХ СЛОЁВ ГАЗА В КОНФИГУРАЦИИ Р. МИЗЕСА

Н.С. Новаковский

аспирант, e-mail: [email protected]

Уральский государственный университет путей сообщения (УрГУПС)

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

Ключевые слова: сильное сжатие газа, конечно-разностный метод «Ромб», метод отслеживания особенностей.

Введение

Математическое описание процесса безударного сжатия идеального газа до любого наперёд заданного значения плотности, в том числе до бесконечной плотности, представляет интерес в связи с проблемой управляемого термоядерного синтеза [1,2]. Идея безударного сжатия весьма привлекательна для получения сколь угодно больших плотностей при минимальных энергетических затратах благодаря отсутствию в течении ударных волн, т.е. сохранению в процессе сжатия начальной энтропии [3]. О привлекательности этой идеи можно судить по большому числу публикаций на эту тему (подробную библиографию см. в [4]).

В [4,5] разработана математическая теория безударного сильного сжатия идеального газа. В частности, для случая сжатия цилиндрически V =1 и сферически V = 2 симметричных слоев политропного газа с показателем 7 > 1 доказано, что непрерывная состыковка двух течений даёт решение задачи о безударном сильном сжатии до любой наперёд заданной плотности ненулевой массы газа.

В работе [6] предложен алгоритм расчёта безударного сильного сжатия в обратном направлении изменения времени одномерных слоёв первоначально

однородного и покоящегося газа с р0 = 1 до любой наперёд заданной конечной постоянной плотности р* > 1 с последующим восстановлением закона движения поршня, сжимающего слой газа. В работе [7] представлены результаты расчётов вышеупомянутым алгоритмом ряда одномерных задач с различной симметрией при сжатии одномерных слоёв снаружи.

В работе [8] предложен алгоритм расчёта безударного сильного сжатия изнутри одномерных слоёв газа при возрастании времени, основанный на конечно-разностном методе «Ромб» [9]. При этом используется закон движения поршня, сжимающего слой газа, восстановленный при решении задачи в обратном направлении изменения времени алгоритмом, предложенным в [6].

Расчёты в прямом направлении изменения времени по схеме «Ромб» выявили области течения, в которых предложенный алгоритм искажает решение в том числе в случае плоской симметрии (V = 0).

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

1. Математическая постановка одномерной задачи безударного сильного сжатия газовых слоёв при возрастании времени

В момент времени £ = £0 дан одномерный плоско- (V = 0), цилиндрически-(V = 1) или сферическисимметричный (V = 2) слой идеального газа. Ширина слоя — ^о, масса т0, плотность р0 = 1. Будем считать, что нам известны (рис. 1): гр(£), ир(£) — траектория движения и скорость поршня, сжимающего этот слой газа, гс+(£) — траектория звуковой характеристики С+, отделяющей область покоящегося газа, обозначенная на рисунке цифрой 0, от области центрированной волны сжатия, обозначенной 1 (эти параметры течения восстанавливаются алгоритмом из работы [8]).

Необходимо построить решение при возрастании времени в тех частях областей 1 и 2, которые заключены между траекторией поршня (обозначенной на рис. 1 штрихованной линией) и характеристикой С+. В области 0 решение известно — это покоящийся газ р = 1,и = 0.

В предыдущей работе автора [8] для построения решения в прямом направлении изменения времени конечно-разностный метод «Ромб» использовался во всей области от поршня до г*. Это приводит к «размазыванию» слабого разрыва, через который стыкуется область однородного покоя и область обобщённой центрированной волны сжатия (см. рис. 2).

Исходя из вышесказанного, будем строить решение методом «Ромб» только в области между поршнем и характеристикой С+.

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

Рис. 1. Начальные данные для расчёта

Рис. 2. Неточности в описании сжатия слоя методом «Ромб»

динамики в лагранжевых массовых переменных в одномерном случае.

дь = ди дЬ дт' ди = _ др дЬ дт'

дЕ + ьд (ри) = 0 ~д¥ дг '

(1)

где Ь'Т'т'Р'У = р,£,Р — соответственно: время, эйлерова координата, массовая лагранжева координата, плотность, удельный объём, удельная внутренняя

и2

энергия, давление, Е = е + --удельная полная энергия.

Система замыкается уравнением состояния идеального газа р = (7 — 1)рс^Т е = с,0Т' где сь'Т — теплоёмкость при постоянном объёме и температура газа соответственно.

К системе (1) присоединяются граничные и начальные условия в следующем виде

u|r=rl (t) = Up(t),

u\r=rrg (t) = 0, (2)

u(x, 0) = 0,

где ri,rrg — левая (left) и правая (right) границы отрезка, на котором решается система (1). В начальный момент времени t = 0 правая и левая границы совпадают, т.к. поршень, стартуя, порождает возмущение в покоящемся газе, распространяющееся со скоростью звука rl = rp(0) = rrg = rC+ (0). В каждый последующий момент времени t = tn : rl = rp(tn) < rrg = rC+ (tn).

Применив для (1) неявную аппроксимацию и введя шаг по времени т, получим:

v

n+i —v- = / du_\n+1

\dm) '

un+1 _ un__ ( dp\n+1

^ dm) , (3)

+ ^ п+1 = 0.

Систему (3) будем решать, разбивая исходный отрезок [гг,ггд] на N ячеек (пронумеруем узлы полученной сетки г0 = гг, г1, г2,..., гм = ггд). Граничные условия на левой и правой границах запишем так:

иП = «Р(£п),

п (4)

ипм = 0.

Для решения системы (3) используем алгоритм метода «Ромб» [9]. В настоящей работе будет использована модификация метода в части интегрирования уравнения для координат и уравнения энергии. Для повышения точности при интегрировании по времени будем вместо узловых значений давления и скорости с итерации ^ + 1 использовать полусуммы:

и = 1 ■ (и^+1 + ип) , Р = 1 ■ (Р^+1 + рп) ,

где ^ — номер итерации по нелинейности давления.

Ещё раз подчеркнём, что при комбинированном методе решения задачи сильного сжатия одномерного слоя газа будем использовать известную траекторию звуковой характеристики С+, а методом «Ромб» решение на каждом

временном слое строится только в области гр(£п),гс+ (£п)

2. Описание алгоритма комбинированного метода решения задачи сильного сжатия

Далее будет описан шаг комбинированного алгоритма. Пусть расчёт доведён до момента времени £п. Расчёт на п-м временном шаге методом «Ромб» вёлся в

первых N интервалах расчётной сетки. Далее:

1) выбирается шаг по времени тп+1, исходя из общих соображений точности расчёта;

2) используя дискретный набор значений ир(Ь), полученный при решении задачи сильного сжатия в обратном направлении изменения времени [6,8], с помощью линейной интерполяции вычисляется ир (¿п + тп+1). В соответствии с известной траекторией гс+ (¿) получаем гс+ (Ьп + тп+1);

3) добавляется к расчётной области ^ + 1)-й интервал, так что ты+1 = гс+ (¿п + тп+1), значения газодинамических параметров на п-м слое в добавленном интервале соответствуют фоновому течению;

4) делаются вычисления согласно формулам метода «Ромб». Таким образом строим решение на новом временном слое.

Метод «Ромб» является неявным и обладает безусловной устойчивостью. Следовательно, можно выбирать временной шаг достаточно произвольно. Тем не менее из условия точности расчёта временные шаги, начиная со второго, выбираются в соответствии с критерием Куранта ^^ ^ 1. Для старта алгоритма, т.е. для обеспечения достаточного количества узлов расчётной сетки внутри рассчитываемой методом «Ромб» области, первый шаг т1 делается большим. Если при этом для сходимости метода требуется более 10 итераций, счёт прерывается и начинается заново с более мелким шагом.

Рис. 3. Добавление расчётного интервала на правой границе

Опишем более подробно пункт 3 предложенного алгоритма (см. рис. 3). На текущем шаге во всей области г^ < г < г* газодинамические параметры соответствуют фоновому течению. При переходе на следующий временной слой правая граница смещается в соответствии с траекторией С+-характеристики. Величины с п-го слоя во всех интервалах кроме ^ + 1)-го берутся равными расчётным значениям, полученным методом «Ромб». Массовая координата ^ + 1)-го интервала вычисляется по фоновой плотности в зависимости от симметрии задачи шм+1 = р0■ (г^Д — г^+1), где V = 0,1, 2 для плоской, цилиндрической и сферической симметрии соответственно.

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

3. Результаты расчётов

Для проверки работоспособности алгоритма сначала был проведён расчёт плоскосимметричной задачи. Были взяты следующие параметры газа: 7 = 1.4, ш* = 10. Точное решение в области между сжимающим поршнем и характеристикой С+ имеет вид центрированной волны сжатия. На рисунке 4 приведены профили скорости и скорости звука для тестовой задачи. Чёрной сплошной линией здесь обозначено точное решение, синей пунктирной — численное. Можно утверждать, что модифицированный метод хорошо описывает область центрированной волны и область слабого разрыва вплоть до момента времени г = 0.998г.

Рис. 4. Результаты расчёта тестовой задачи на различные моменты времени

Далее было проведено сравнение на тестовой задаче расчётных значений, полученных методом «Ромб», без выделения особенностей, и комбинированным методом. Из рисунка 5 можно сделать следующий вывод: представляемый метод описывает решение в области волны сжатия и в области слабого разрыва принципиально точнее исходного метода «Ромб».

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

г

Рис. 5. Сравнение расчётов двумя методами

баланса по следующей формуле:

Е шк • (ЕП - + ж

дЕп = ^-, (5)

Е Шк • Ек

к= 1'М

где т,ЕП,Е0 — масса, удельная полная энергия интервала на текущий и на начальный момент времени, Ж — работа сжимающего поршня, вычисляемая по формуле:

г"

ж = ^ 1(р1-1 + р1 )(у1 - V1-1) «у Р ¿V,

1=2 п г0

где Р1,У1 — давление на левой границе и объём системы на 1-м временном шаге, полученные по вычисленным ранее значениям скорости звука на поршне и координаты поршня соответственно.

Дисбаланс энергии на различные моменты времени для одного из расчётов приведён в таблице 1.

Таблица 1. Дисбаланс энергии на различные моменты времени

0.95** 0.975** 0.99** 0.995** 0.998**

-5.0861 • 10-5 -8.7531 • 10-5 -2.9584 • 10-4 -3.3866 • 10-4 -4.6709 • 10-4

Дисбаланс связан с неточностью определения скорости поршня на текущем шаге по времени с помощью линейной интерполяции. Из таблицы 1 видно, что он практически не меняется с начала расчёта и его величина незначительна.

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

В качестве эталонного используется решение соответствующих задач в обратном направлении изменения времени методом характеристик с пересчётом.

Рис. 6. Результаты расчёта на различные моменты времени 7 = 1.4, ш* = 10, V =1

Стоит отметить, что область стыковки волны сжатия и покоящегося газа описывается достаточно точно. Для примера приведены два решения различных задач. На рисунке 6 приведены профили искомых величин на различные моменты времени для цилиндрическисимметричной задачи 7 = 1.4, ш* = 10, V =1, а на рисунке 7 — для сферическисимметричной задачи 7 = 5/3, ш* = 10, V = 2.

Рис. 7. Результаты расчёта на различные моменты времени 7 = 5/3, ш* = 10, V = 2

4. Выводы

Представленные результаты позволяют сделать некоторые выводы.

1. Расчёты при возрастании времени неявным конечно-разностным методом «Ромб» демонстрируют работоспособность предлагаемой методики. При этом область волны сжатия адекватно описывается при приближении к моменту времени t*. Точность расчётов не зависит от симметрии задачи. Основная погрешность расчёта локализована в достаточно узкой области вокруг слабого разрыва, отделяющего волну сжатия от покоящегося газа.

2. Реализованная модификация метода «Ромб» обеспечивает отслеживание слабого разрыва. Точность расчёта становится принципиально лучше.

3. При приближении к моменту времени t* в области слабого разрыва наблюдается опережение нижним краем фронта волны сжатия звуковой бихарактеристики, что не позволяет проводить расчёт далее. При этом градиент искомых функций достигает значения 45.

4. Для получения в расчётах при возрастании времени значительных масс сжатого газа требуется привлечение дополнительных аналитических методов.

Автор выражает своему научному руководителю профессору С.П. Баутину

признательность за внимание, помощь и поддержку.

Литература

1. Забабахин Е.И., Забабахин И.Е. Явления неограниченной кумуляции. М. : Наука, 1988. 173 с.

2. Накколс Дж.Г. Осуществимость инерциально-термоядерного синтеза // Успехи физ. наук. 1984. Т. 143, № 3. С. 467-482.

3. Долголева Г.В., Забродин А.В. Кумуляция энергии в слоистых системах и реализация безударного сжатия. М. : ФИЗМАТЛИТ, 2004. 71 с.

4. Баутин С.П. Математическое моделирование сильного сжатия газа. Новосибирск : Наука, 2007. 308 с.

5. Баутин С.П. Математическая теория безударного сильного сжатия идеального газа. Новосибирск : Наука, 1997. 160 с.

6. Баутин С.П., Николаев Ю.В. Об одном методе расчёта безударного сильного сжатия одномерных слоёв газа // Вычислительные технологии. 2000. Т. 5, № 4. С. 312.

7. Николаев Ю.В. О численном решении задачи безударного сильного сжатия одномерных слоёв газа // Вычислительные технологии. 2001. Т. 6, № 2. С. 104-108.

8. Новаковский Н.С. Математическое моделирование сильного сжатия одномерных слоёв газа в конфигурации Р. Мизеса // Математические структуры и моделирование. 2016. № 3(39). С. 93-109.

9. Гаджиев А.Д., Писарев В.Н. Неявный конечно-разностный метод «Ромб» для численного решения уравнений газовой динамики с теплопроводностью. // Ж. выч. ма-тем. и матем. физ. 1979. Т. 19, № 5. C. 1288-1303.

THE COMBINED NUMERICAL METHOD FOR SOLVING THE ONE-DIMENSIONAL IDEAL GAS SHOCK-FREE STRONG COMPRESSION PROBLEM IN R. MISES CONFIGURATION

N.S. Novakovskiy

Graduate Student, e-mail: [email protected]

Ural State University of Railway Transport (USURT)

Abstract. A method for solving the one-dimensional ideal gas shock-free strong compression problem in R. Mises configuration is proposed. The method combines finite-difference method "ROMB" and tracking feature method. The method allows to calculate gas-dynamic characteristic (velocity, density, etc.) of ideal gas layer while time increases and provides better accuracy in comparison with other finite-difference methods. The accuracy of the proposed method was demonstrated in calculations of test plane-symmetry problem. Exact solution and numerical one agree quite well. Numerical results of solving one-dimensional problems with different symmetry and gas characteristic are also shown. The main results of numerical simulations are shown in graphs and tables.

Keywords: gas strong compression, finite difference method "Romb", discontinuity tracking method.

Дата поступления в редакцию: 27.12.2016