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

УДК 519.63+532.593

DOI: 10.14529/ mmp190205

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

Ю.М. Ковалев1, П.А. Кузнецов1

1 Южно-Уральский государственный университет, г. Челябинск, Российская Федерация

В данной работе предложена модификация метода крупных частиц. Проведен численный анализ различных модификаций метода крупных частиц применительно к задачам волновой динамики (газовой динамики). Решены задачи расчета распада произвольного разрыва, а также задачи о распространении стационарных ударных волн с отражением от жесткой стенки. Было показано, что предложенная модификация метода крупных частиц наилучшим образом совпадает с аналитическим решением задачи об отражении плоской ударной волны от жесткой стенки. Проведенный численный анализ показал, что данная модификация позволяет проводить устойчивые расчеты течений с большими градиентами изменения параметров. Значительным достоинством предложенной модификации является тот факт, что рассмотренные в работе задачи могут быть решены без введения в законы сохранения «искусственной» вязкости, а также при больших числах Куранта.

Ключевые слова: ударная волна; распад произвольного разрыва; число Куранта; искусственная вязкость; крупные частицы.

Введение

В настоящее время математическое моделирование стало полноценным инструментом при проведении исследований широкого круга задач газовой динамики, а при изучении быстропротекающих процессов математическое моделирование очень часто является единственным способом получения углубленной и достоверной информации о многих физических и химических особенностях исследуемых течений газа. Для проведения численных экспериментов необходимо разработать и адаптировать необходимый вычислительный аппарат, который бы позволил решать широкий круг практических задач. Одним из наиболее удобных и универсальных численных методов является метод крупных частиц [1] (МКЧ), использующий принцип расщепления по физическим процессам [2-4]. Метод крупных частиц имеет достаточно простой алгоритм, что определило его широкое распространение. Однако при математическом моделировании течений с большими градиентами изменения параметров приходится вводить искусственную вязкость [5-7] и снижать сеточное число Куранта.

Целью настоящего исследования является анализ возможностей модификаций МКЧ, приведенных в работах [5,8], с модификацией, предложенной авторами данной статьи, при решении задач, представляющих интерес для изучения быстропротека-ющих процессов:

1) распространение ударных волн;

2) отражение стационарной ударной волны от жесткой стенки;

3) расчет распада произвольного разрыва.

1. Постановка задачи и описание метода

Рассмотрим одномерное движение идеального сжимаемого газа, описываемое дифференциальными уравнениями Эйлера в дивергентном виде (уравнения неразрывности, импульса и энергии)

др

— + сиу(ри) = О,

дри др

¥ + (11¥(Н + ^ = 0, (1)

др е

—Ь с11у(рЕи) + сИу(рм) = 0.

Замыкается система (1) уравнением состояния и дополняется уравнением удельной полной энергии газа

Р = ер(7 - 1), и2

Здесь р, р, и, Е, е, Ь, х, 7 - давление, плотность, скорость, удельная полная и внутренняя энергии, время, декартова координата и показатель адиабаты Пуассона соответственно. В соответствии с решаемыми задачами к системе законов сохранения добавляются начальные и граничные условия.

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

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

II - лагранжев этап, где при движении жидкости вычисляются потоки массы через границы эйлеровых ячеек;

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

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

Ах рЧ у '

Ёп = Еп - ~ -^-1/2^-1/2

* * Дж р™' [ }

где р, р, и, Е, Ь, х - соответственно давление, плотность, скорость, полная энергия, время и декартова координата.

Величины с дробными индексами, относящиеся к границам ячеек, находятся так:

/i in I /11n <rtn I <rin

ra _ Uí T »¿+1 ra _ Pi ^ Pí+1 Ui+1/2 — «л > Pi+1/2 ~

2 ' 1г+1/2 2

В формулах (2) и (3) вместо р можно использовать р + д. Здесь д - псевдовязкость.

В работе [5] автор предлагает введение «локальной> псевдовязкости, для минимизации недостатков метода. Она представляет собой комбинацию линейной и квадратичной псевдовязкости, значение которой определяется на границах ячеек после анализа сравнения абсолютных значений первой и второй производных решения:

ди дх

г л du

q = —орАх—

CTia + a2Ax

5 = const > 0, 0 < a1 < 1,

0, ¿1,

¿2,

д2 u

дх2 д2 и

Аж - g IIм! < 0

0 < a2 < 1,

ди

дх2 д2и

Ax — g

дх

IttI < 0

дх

л > 0, дх

ди

дх2

дх

Я < о,

дх

> 0, g > 0.

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

(4)

-5рпг+1 «+1 - <)(<7ia^+i + <j2|<+1 - <|),

где

пп

{Рг+1 - Р?) , CI = Ö (РШ - Рг)

2<¿

„га Pí+h

0, |u n- 2ui+1 + Ui+2 1 |< g| Un - - Un | ui+2 1 , un+1 — un > 0,

5 = < 5b — 2ui+1 + Ui+2 | < g |un+1 — Un | ui+2| , un+1 — un < 0,

( 52, 0„n 1 ».п. | 2ui+1 + ui+2| >g |Un 1 ui+1 — Un | g = const > 0.

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

Отметим, что недостатком введения псевдовязкости на эйлеровом этапе является необходимость определения констант ^2, 9 в формуле (4) для каждой

конкретной задачи.

5

1

2

2. Модификации метода крупных частиц

В работе [10] был проведен сравнительный анализ различных модификаций МКЧ [8,9]. При решении тестовых задач лучшие результаты были получены при использовании модификации МКЧ, предложенной Ю.А. Гришиным и В.А. Зенки-ным [8]. В данной модификации предлагаются иные формулы для расчета предварительного давления, скорости и полной энергии на эйлеровом этапе.

Предварительные значения давления на n +1 временном слое для i — 1/2 и i + 1/2 границ между ячейками i — 1, i и i + 1 определяются по формуле:

Plíi/2 = ^V^ (l - (7 - 1)«+1 - • (5)

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

К+1 = <-Рг+1/\ Рг~1/2 —• (6)

Ах рЧ у '

После этого вычисляются скорости на границах ячеек по формуле

Щ + иг+1

и

4+1/2 2

Формула (3) при этом заменяется на следующую:

Рт+1 = ™ _ - Р7-1/2и7-1/2 Аг

Ах р?'

Для повышения устойчивости счета МКЧ в зонах нулевой скорости авторами данной статьи предлагается следующая модификация метода. На эйлеровом этапе для определения давления на границах ячеек в зонах сжатия вместо формулы (5) следует использовать соотношения на ударной волне, по аналогии с методом В.Ф. Куропа-тенко [11]. Если и™+1 — и™ < 0, то

р = Г РП+1 + — <+1) , р™ <Р1+1, (8)

Рг+-2 \ р?-\¥г(и?-и?+1), &>р?+1, [ ]

где Wi - скорость сеточного разрыва.

После определения предварительных скоростей в серединах ячеек по формуле (6) скорости на границах ячеек в зонах сжатия считаются заданными скоростями за разрывом. Если и™+1 — и™ < 0, то

и 1 = { <> Р?<Р?+ г,

Предварительные значения полной энергии в серединах ячеек определяются по формуле (7). Лагранжев и заключительный этапы проводятся аналогично базовому МКЧ.

3. Распад произвольного разрыва

Рассмотрим задачу распада произвольного разрыва (РПР). Бесконечная труба, заполненная воздухом (7 = 1,4), в точке х = 0 разделена перегородкой. При Ь = 0 состояние газа описывается следующей функцией:

= Г (рь 0 рь), х < 0, \ (рп 0 рп), х > 0.

г=о

Здесь рь, рь, рп, рп - плотность и давление до перегородки и за перегородкой соответственно.

Пусть

1 0 1

0,125 0 0, 1

(9)

В момент времени Ь = 0 перегородка мгновенно убирается, и эта конфигурация распадается на волну разрежения, распространяющуюся влево, и ударную волну с контактным разрывом, распространяющиюся вправо.

На рис. 1 представлены графики численного решения задачи (9) предлагаемой авторами модификацией МКЧ, модификацией МКЧ [8], МКЧ с использованием псевдовязкости, а также график аналитического решения, полученного из соотношений для распада произвольного разрыва, на момент времени Ь = 0, 502, с сеточным числом Куранта К = 0,1. Количество ячеек N = 100. На данных рисунках пунктирная линия соответствует модификация МКЧ [8], штрих-пунктирная - базовый МКЧ с псевдовязкостью (4), сплошная линия - модификация МКЧ, предлагаемая авторами данной работы, штриховая линия - аналитическое решение. На рис. 1 видно, что полученные численные решения в целом совпадают с аналитическим решением.

Рис. 1. Решение задачи РПР. Графики давления (слева) и скоростей (справа)

4. Распространение стационарной ударной волны

Рассмотрим распространение стационарной ударной волны (число Маха М = 2, 3). Параметры среды в начальный момент времени Ь = 0: р =1, и = 0, Е = 0 -в области невозмущенного течения, а на границе определялись из соотношений для ударной волны.

На рис. 2 представлены графики стационарных ударных волн (УВ), полученные при помощи предлагаемой авторами модификацией МКЧ, модификацией МКЧ [8], МКЧ с использованием псевдовязкости, а также график аналитического решения на момент времени Ь = 0, 2361 с числом Куранта К = 0,1. Количество ячеек N = 100. Обозначения аналогичны рис. 1.

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

Рис. 2. Профили давлений (слева) и скростей (справа) стационарной УВ

5/"~ч о о о

. Отражение стационарном ударной волны от жесткой стенки

Рассматривается та же задача, что и в разделе 4. Ударная волна доходит до жесткой стенки (ж = 1) и отражается от нее.

На рис. 3 представлены графики расчетных отраженных стационарных ударных волн (УВ), а также график аналитического решения на момент времени £ = 0, 3879 с числом Куранта К = 0,1. Количество ячеек N = 100. Обозначения аналогичны рис. 1.

0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 х 0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 х

Рис. 3. Профили давлений (слева) и скоростей (справа) отраженной стационарной УВ

Из рис. 3 видно, что при решении задачи модифицированным МКЧ [8] возникают сильные осцилляции за фронтом УВ. При использовании псевдовязкости (4) осцилляции гасятся. Однако, как уже отмечалось выше, использование псевдовязкости требует подбора эмпирических констант для каждой конкретной решаемой задачи. На этом фоне более универсальной является модификация с расчетом предварительных значений эйлерового этапа на ударной волне из соотношений Гюгонио. В этом случае нет необходимости подбирать эмпирические константы. Стоит отметить, что ударный фронт оказывается на 2-3 ячейки шире по сравнению с другими проведенными расчетами. Добавим, что предложенная авторами модификация позволяет получать устойчивые расчеты при больших числах Куранта - до К = 2 для решаемых в данной статье задач. Это большим преимуществом модификации при решении более сложных задач, требующих для расчета больших временных затрат.

Заключение

Проведенный в данной работе численный анализ различных модификаций [5, 8] метода крупных частиц показал, что предложенная выше модификация эйлерова этапа, использующая алгоритм, представленный в работе [11], позволяет проводить устойчивые расчеты течений с большими градиентами изменения параметров с числом Куранта равным 2 без введения псевдовязкости, хорошо совпадающие с аналитическими решениями.

Литература

1. Белоцерковский, О.М. Метод крупных частиц в газовой динамике / О.М. Белоцерков-ский, Ю.М. Давыдов. - М.: Наука, 1982.

2. Ковеня, В.М. Алгоритмы расщепления при решении многомерных задач аэродинамики / В.М. Ковеня. - Новосибирск: Изд-во СО РАН, 2014.

3. Грищенко, Д.С. Модификация метода крупных частиц для исследования течений газовзвесей / Д.С. Грищенко, Ю.М. Ковалев, Е.А. Ковалева // Вестник ЮУрГУ. Серия: Математическое моделирование и программирование. - 2015. - Т. 8, № 2. - С. 36-42.

4. Ковалев, Ю.М. Анализ некоторых модификаций метода крупных частиц на примере исследования течений газовзвесей / Ю.М. Ковалев, Е.А. Ковалева, Е.Е. Пигасов // Вестник ЮУрГУ. Серия: Математика. Механика. Физика. - 2015. - № 3. - С. 71-77.

5. Ивандаев, А.И. Об одном способе введения «псевдовязкости» и его применении к уточнению разностных решений уравнений гидродинамики /А.И. Ивандаев // Журнал вычислительной математики и математической физики. - 1975. - Т. 15, № 2. - С. 523-527.

6. Губайдуллин, А.А. Применение модифицированного метода «крупных частиц» к решению задач волновой динамики / А.А. Губайдуллин, А.И. Ивандаев // Журнал вычислительной математики и математической физики. - 1976. - Т. 16, № 4. - С. 1017-1026.

7. Губайдуллин, А.А. Модифицированный метод «крупных частиц» для расчета нестационарных волновых процессов в многофазных дисперсных средах / А.А. Губайдуллин, А.И. Ивандаев, Р.И. Нигматулин // Журнал вычислительной математики и математической физики. - 1977. - Т. 17, № 6. - С. 1531-1544.

8. Гришин, Ю.А. Повышение устойчивости вычислительного алгоритма метода крупных частиц / Ю.А. Гришин, В.А. Зенкин // Наука и образование. - 2011. - № 13. - С. 41-45. - URL: http://old.technomag.edu.ru/doc/221488.html

9. Гришин, Ю.А. Новые расчетные схемы на базе метода крупных частиц для моделирования газодинамических задач / Ю.А. Гришин, В.Н. Бакулин // Доклады Академии наук. - 2015. - Т. 465, № 5. - С. 545-548.

10. Kovalev Yu.M. Analysis of Some Modifications of the Large-Particle Method to Model Wave Dynamics Problems / Yu.M. Kovalev, P.A. Kuznetsov // Journal of Computational and Engineering Mathematics. - 2018. - V. 5, № 3. - P. 38-48.

11. Куропатенко, В.Ф. Метод расчета ударных волн / В.Ф. Куропатенко // Доклады Академии наук СССР. - 1960. - Т. 133, № 4. - C. 771-773.

Юрий Михайлович Ковалев, доктор физико-математических наук, профессор, заведующий кафедрой «Вычислительная механика:», Южно-Уральский государственный университет (г. Челябинск, Российская Федерация), kovalevym@susu.ru.

Павел Анатольевич Кузнецов, аспирант, кафедра «Вычислительная механика:», Южно-Уральский государственный университет (г. Челябинск, Российская Федерация), kuznetcovpa@susu.ru.

Поступила в редакцию 8 ноября 2018 г.

MSC 76L04, 76M20 DOI: 10.14529/mmp190205

MODIFICATION OF THE LARGE-PARTICLE METHOD TO SOLVE SHOCK WAVES AND RAREFACTION WAVES PROPAGATION

Yu.M. Kovalev1, P.A. Kuznetsov1

1 South Ural State University, Chelyabinsk, Russian Federation

E-mails: kovalevym@susu.ru, kuznetcovpa@susu.ru

In this article, we present a modification of the large-particle method. We perform a numerical analysis of various modifications of the large-particle method applied to problems of wave dynamics (gas dynamics). We solve the problems on calculation of the decay of an arbitrary discontinuity, as well as the problems on propagation of stationary shock waves with reflection from a solid wall. Calculations show that the solution obtained by the modified large-particle method best coincides with the analytical solution to the shock wave reflection problem from a solid wall. The numerical analysis shows that this modification allows to carry out stable flow calculations with large parameter gradients. A significant advantage of this modification is the fact that the considered problems can be solved without introducing the "artificial" viscosity into the laws of conservation, and with large Courant numbers.

Keywords: shock wave; decay of an arbitrary discontinuity; Courant number; artificial viscosity; large particles.

References

1. Belocerkovskii O.M., Davydov Ju.M. Metod krupnyh chastic v gazovoy dinamike [Large-Particle Method in Gas Dynamics]. Moscow, Nauka, 1982. (in Russian)

2. Kovenia V.M. Algoritmy rasshhepleniya pri reshenii mnogomernyh zadach ayerodinamiki [Numerical Method of Particles in Cells for Hydrodynamics Problems: Computing Methods in Hydrodynamics]. Novosibirsk, Publishing House of Siberian Branch of the Russian Academy of Sciences, 2014. (in Russian)

3. Grishchenko D.S., Kovalev Yu.M., Kovaleva E.A. Modification of Method of Large Particles for Research of Currents of Gas-Suspensions. Bulletin of the South Ural State University. Series: Mathematical Modelling, Programming and Computer Software, 2015, vol. 8, no. 2, pp. 36-42. DOI: 10.14529/mmp150203 (in Russian)

4. Kovalev Yu.M., Kovaleva E.A., Pigasov E.E. The Analysis of Some Modifications of the Large-Particles Method on the Basis of Research of Gas-Suspension Currents. Bulletin of the South Ural State University. Series: Mathematics. Mechanics. Physics, 2015, vol. 7, no. 3, pp. 71-77. (in Russian)

5. Ivandaev A.I. A Method of Introducing "Pseudoviscosity" and Its Use for Improving the Difference Solutions of Hydrodynamic Equations. USSR Computational Mathematics and Mathematical Physics, 1975, vol. 15, no. 2, pp. 238-242. DOI: 10.1016/0041-5553(75)90063-4

6. Gubaidullin A.A., Ivandaev A.I. Use of a Modified Large-Particle Method for Solving Problems of Wave Dynamics. USSR Computational Mathematics and Mathematical Physics, 1976, vol. 16, no. 4, pp. 196-205. DOI: 10.1016/0041-5553(76)90018-5

7. Gubaidullin A.A., Ivandaev A.I., Nigmatulin R.I. A Modified "Coarse Particle" Method for Calculating Non-Stationary Wave Processes in Multiphase Dispersive Media. USSR Computational Mathematics and Mathematical Physics, 1977, vol. 17, no. 6, pp. 180-192. DOI: 10.1016/0041-5553(77)90183-5

8. Grishin Yu.A., Zenkin V.A. Improving the Stability of the Numerical Algorithm of the Method of Large Particles. Journal of Science and Education, available at: www.old.technomag.edu.ru/doc/221488. (in Russian)

9. Grishin Yu.A., Bakulin V.N. New Calculation Schemes Based on the Large-Particle Method for Modeling Gas-Dynamic Problems. Doklady Physics, 2015, vol. 60, no. 12, pp. 555-558. DOI: 10.7868/S0869565215350091

10. Kovalev Yu.M., Kuznetsov P.A. Analysis of Some Modifications of the Large-Particle Method to Model Wave Dynamics Problems. Journal of Computational and Engineering Mathematics, 2018, vol. 5, no. 3, pp. 38-48. DOI: 10.14529/jcem180304

11. Kuropatenko V.F. [A Shock Calculation Method]. Doklady Physics, 1960, vol. 3, no. 4, pp. 771-772. (in Russian)

Received November 8, 2018