CC BY
105
www.volsu.ru
ПРИКЛАДНАЯ МАТЕМАТИКА
DOI: http://dx.doi.oгg/10.15688/j•volsu1.2015.1.2
УДК 519.6 ББК 22.193
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ЗВУКОВЫХ ВОЛН В АКТИВНЫХ СРЕДАХ С ИСПОЛЬЗОВАНИЕМ СХЕМЫ MUSCL1
Бочкарева Екатерина Викторовна
Студентка Института математики и информационных технологий, Волгоградский государственный университет [email protected], [email protected]
просп. Университетский, 100, 400062 г. Волгоград, Российская Федерация
Храпов Сергей Сергеевич
Кандидат физико-математических наук, доцент кафедры информационных систем и компьютерного моделирования, Волгоградский государственный университет [email protected], [email protected]
просп. Университетский, 100, 400062 г. Волгоград, Российская Федерация
Аннотация. На основе реализации метода ИЬЬС и схемы МиБСЬ исследована неустойчивость акустических волн в активных средах. Проведено сравнение метода типа Годунова с использованием численных схем первого и второго порядка точности. Исследовано влияние ограничителей наклонов и ю способов вычисления потоков на качество численного решения. § Ключевые слова: численные схемы, ограничители наклонов, метод типа
Годунова, МШСЬ, ИЬЬС.
и
§ Введение
с га а
^ В данной работе подробно рассмотрен метод ИЬЬС, проведено сравнение точно-оз сти различных приближенных методов решения задачи Римана о распаде произволь-^ ного разрыва. Решение этой задачи необходимо для построения численных методов, £ принадлежащих типу Годунова. Для пересчета величин по времени использовалась схе-р ма типа МиБСЬ, которая предполагает простейший двушаговый пересчет, называемый предиктор-корректор. На основе этого алгоритма была смоделирована динамика волно-@ вых процессов в активных средах.
1. Активная среда
В реальных средах звуковые волны затухают вследствие вязкости среды и молекулярного рассеяния. Однако в активных средах распространение акустических волн происходит иначе.
Газодинамика термодинамически неравновесных сред, таких как колебательно-возбужденный газ, неизотермическая плазма, химически активные смеси, среды с неравновесными фазовыми переходами и т. д., является в современном мире предметом экспериментальных исследований и теоретических работ. В таких средах коэффициент второй (объемной вязкости) и акустическая дисперсия могут быть отрицательными. Акустическая активность сред приводит к сильному изменению структуры и динамики нелинейных волн, в ней распространяющихся. Многочисленные эксперименты, проведенные в неравновесных средах, свидетельствуют о существенном изменении структуры газодинамического возмущения. А именно, в неравновесных газоплазменных средах наблюдается расщепление фронта ударной волны, изменение ее скорости, амплитуды и ширины фронта, образование предвестников. Кроме того, низкочастотный коэффициент гидродинамической нелинейности в таких средах является сложной функцией стационарной степени неравновесности и может быть отрицательным.
Эти новые акустические свойства необходимо учитывать при рассмотрении различных газодинамических явлений, протекающих в этих средах. Прежде всего их влияние может быть заметным при исследовании структуры газодинамических возмущений малой, но конечной амплитуды [6].
Таким образом, активная среда — вещество, в котором распределение частиц (атомов, молекул, ионов) по энергетическим состояниям не является равновесным и хотя бы для одной пары уровней энергии осуществляется инверсность населенностей. То есть, для того чтобы среда была неравновесной, необходимо избирательное возбуждение атомов или молекул, обеспечивающее избыточное заселение одного или нескольких верхних уровней энергии по сравнению с нижележащим уровнем.
2. Численная схема
2.1. Уравнения газодинамики в активной среде
Математически описать распростанение звуковых волн в среде можно с помощью уравнений Эйлера. Если среда является неравновестной, то система уравнений Эйлера в одномерном случае имеет вид:
др д т + дX(р и)
0,
дри д . 2. др . .
~дГ + д:X {Ри ) = - дХ +р1(Х,1),
дЕ д д
~Ж + дХ(иЕ) = - дХ(ир)+р [^(р,р) - 1(р,(x,
дЕи +Х ( Е ) *
"Ж" + дХ(иЕ ) = Р
(1)
Ее Е + Я(р, р)
т(p, р)
где Я(р,р), 1(р,р) — мощности источника накачки и теплоотвода (<5 = I = з/т0, 0.5 <
< s < 1.5); Ev — колебательная энергия; Ее — равновесное значение колебательной энергии; т(р,р) — время колебательной релаксации. Для активной среды энергия находится по формуле
2
Е = + + Е„, (2)
7 — 1 2
где 7 — показатель адиабаты.
Колебательная энергия для неравновесной среды определяется через ее равновесное значение, время колебательной релаксации и мощность источника накачки:
Еи = Ее + т('ро, po)Q(po, ро), (3)
где ро, Ро — стационарные значения плотности и давления.
Время колебательной релаксации вычисляется по формуле [4; 5]
r(p. Р) = М JL Г-(2. ),т, То = С. = Ж (4)
Ро Ро 2пCs V Р
где тг — время вращательной релаксации; тт — время поступательной релаксации.
2.2. Метод типа Годунова
Для построения численной схемы воспользуемся явным методом типа Годунова [12]
иГ1 - Ц? + Е+1/2 - Ег-1/2 = 0 (5)
Аг Ах '
где
{р\ / ри
Ц = I ри I ' Е = I ри2 +р | . (6)
\и(Е +р)/
Значения потоков 1 на границах ячеек вычисляются из решения соответствующей задачи Римана. Будем вычислять эти значения различными способами: методом Лакса — Фридрихса (ЬЕ), Хартена — Лакса — ван Лира (ИЬЬ) и модифицированным методом Хартена — Лакса — ван Лира (ИЬЬС).
2.3. Метод HLLC
Метод ИЬЬС предполагает расчет 1 исходя из системы [11]:
Еь, если 0 < БЬ'
Е*ь' если Бь < 0 < 5*' Е^ 1 = (7)
Е*Я' если < 0 < 3Я'
Ед' если 0 > Бя'
где Е*д и Е*^ (индексы К и Ь заменяет К) вычисляются следующим образом:
Б*(8к Ик — Ек) + вкРькБ*
Е
*к
вк — 5*
Рьк =^[рь + Рк + Рь(Бь — иь)(Б* — и£) + Рк(вк — ик )(в* — ия)], (8)
Б* = [0,1,5*].
Величина 5* определяется по формуле
Рк — Рь + Рьиь(вь — иь) — ркик(8к — ик)
5
*
Рь( вь — иь) — рк(3и — ип (9)
вь = иь — аьдь, = ик + акдк,
где
{1, если р* < рк,
1 + ^ (Р */р к — 1)
, если р* > рк.
Аккустическое приближение давления * находится как
р* = тах(0, ^(рь + ря) — ^(иь — ия)ра), Р=2>( Рь + Рк), а=^(аь + ая).
(10)
(11)
2.4. Схема MUSCL
Построение схемы высокого порядка точности осуществляется путем сочетания использования кусочно-линейной аппроксимации величин внутри ячеек с различными алгоритмами пересчета по времени. В данной работе будем использовать схему типа MUSCL, в основе которой лежит простейший двухшаговый пересчет, называемый предиктор-корректор.
Предиктор: первый шаг. Предположим, что внутри дискретных ячеек для всех значений сеточных функций заданы кусочно-линейные распределения вида:
И(Г,хг) = и + Q?(X — Хг), X е X — 1Лх,хг + ±Лх], (12)
где Хг — пространственная координата центра ячейки с номером г, а Q™ — вектор наклонов распределения функции и внутри ячейки в направлении оси х.
Запишем разностное уравнение для учета изменения И по времени в виде:
и г1 — И* + е(и + \Лхдр — е(и — \Лхд?) =0 (13)
лг Лх .
Шаг по времени At определяется из условия устойчивости численной схемы [1; 2]:
. / Ах
At = K min (.—,-—),
¿=1,^ |<| + С/'
(14)
/
а = 4/ У*
■ Р?
где К — критерий Куранта, 0 < К < 1.
Предиктор: второй шаг. Вычислим значение функции и на промежуточном слое по времени Ь + 1АЬ:
и™+1/2 = 2(ИП+1 - ИТ). (15)
Корректор. На шаге «корректор» применим конечно-объемную схему
ИТ+1 - И + Г-+1/2 - ^ =0 (16)
А* Ах
где все значения Р^+1/2 определяются решением задачи Римана со следующими кусочно-постоянными начальными данными:
Ги™+1/2 + 1 АхО™, при хг+1/2 < 0, \и™+1/2 - 2АхОТ+1, при жг+1/2 > 0. Значения потоков Р^-1/2 определяются аналогично.
2.5. Ограничители наклонов
Одним из простых способов вычисления наклонов О™ в дискретной ячейке с номером г для сеточной функции И, где г = 1, 2,..., который обеспечивает устойчивость схемы, является использование функции штшс^ [3; 8]. Наклоны в ячейках рассчитываются по формуле:
И - И-1
г =-,
Иг+1 - И/ (18)
minmсd(r) = тах[0, тт(1, г)].
Функция minmсd выбирает наклон с минимальным значением абсолютной величины, при условии что знаки обоих аргументов совпадают. При различных знаках аргументов функция равна нулю.
Однако существуют и другие ограничители наклонов: — ограничитель ван Лира [9]:
г + |г|
1^16^(7-) = ; (19)
1 + |г|
— ограничитель ван Альбады:
у2 | у
limiter ^ а(г) = -—-; (20)
г2 + 1
— ограничитель «superbee» [10]:
superbee(r) = max[0,min(2г, 1), min(г, 2)]. (21)
3. Проведение вычислительных экспериментов 3.1. Тестирование
В качестве тестовой рассмотрим задачу о распаде произвольного разрыва газа с показателем адиабаты 7 = 1, 4 в точке х0 = 0 с начальными условиями и = 0, р = 1, р = 1, 0 при х < х0, иначе р = 0,1. Определим критерий Куранта К = 0,5. Вычисления проводились на сетках с различным количеством ячеек N = 100; 200; 400; 800. На рисунке 1 приведены результаты тестирования при N = 800.
Рис. 1. Результаты тестирования с использованием метода ИЬЬС, схема первого (а) и второго (б) порядка точности в момент времени £ = 0, 5. Сплошной линией показано аналитическое решение для профиля плотности, ромбами — численное
Таблица 1
Погрешность вычислений при различных методах
Первый порядок точности Второй порядок точности
N ЬР ИЬЬ ИЬЬС ЬР ИЬЬ ИЬЬС
100, Т = 0.5 0.12124 0.099639 0.089542 0.058916 0.048838 0.046978
200, Т = 0.25 0.06090 0.050060 0.044980 0.029604 0.024540 0.023605
400, Т = 0.125 0.030528 0.025092 0.022546 0.014839 0.012300 0.011832
800, Т = 0.0625 0.015283 0.012562 0.011287 0.007428 0.006158 0.005923
Таблица 2
Погрешность вычислений при различных функциях-ограничителях
N штшоб ван Лир ван Альбада эирегЬее
200, Т = 0.5 0.02940 0.022804 0.025170 0.017009
800, Т = 0.5 0.012072 0.007685 0.008751 0.005682
Для каждого из вариантов схемы была рассчитана относительная ошибка по формуле
1 М
6 = ^Е ^ - ^ (22)
г=1
где ра — аналитическое решение; рп — численное. Результаты представлены в таблицах 1, 2.
3.2. Газодинамическое моделирование в активной среде
Для моделирования динамики акустических волн в неравновесной среде зададим вынужденные колебания от источника с начальными условиями и = 0, р = 1, р = = 1, критерий Куранта К = 0, 5. Для вычисления времени колебательной релаксации определим значения тг = -1,0 и тт = -1, 5. Длина волны А = 0,1.
Рассмотрим линейную стадию распространения звуковых волн, для этого зададим Б = 0, 5 и сравним полученный результат с функцией
/ = ^е х '
(23)
где А — амплитуда начального возмущения. На рисунке 2 представлены результаты при N = 800, для функции (23) задан параметр а = 0,43. На линейной стадии хорошо прослеживается экспоненциальный рост амплитуды возмущения.
Для моделирования нелинейной стадии звуковых волн был задан параметр Б = = 1, а = 1, 034. На рисунках 3, 4 представлены распределения возмущения скорости, давления и плотности в звуковой волне. В момент времени £ = 0, 8 начальные импульсы эволюционируют в возмущения типа «ступенька» [4]. Амплитуда, скорость и форма таких возмущений определяются только параметрами неравновесной среды.
Рис. 2. Распределение возмущения скорости в акустической волне на линейной стадии в момент времени £ = 0, 8. Сплошной линией показана огибающая возмущения (23), линией с ромбами — численное решение
Рис. 3. Распределение возмущения скорости в акустической волне на нелинейной стадии в момент времени £ = 0, 8. Сплошной линией показана огибающая возмущения (23), линией с ромбами — численное решение
Заключение
Полученные результаты свидетельствуют о том, что применение схемы MUSCL дает уменьшение погрешности вычислений в 2 раза в независимости от метода решения задачи Римана. Среди ограничителей наклонов наиболее точным является supeгbee. На основе реализации метода HLLC и схемы MUSCL была исследованна неустойчивость акустических волн в активных средах. Полученные результаты свидетельствуют о существенном влиянии стационарной неравновесности (через формируемые ею новые вязкостно-дисперсионные и нелинейные свойства среды) на эволюцию газодинамического возмущения.
Рис. 4. Распределение возмущения давления (а) и плотности (б) в акустической волне на нелинейной стадии в момент времени £ = 0, 8
ПРИМЕЧАНИЕ
1 Работа выполнена при поддержке гранта РФФИ № 15-02-06204.
СПИСОК ЛИТЕРАТУРЫ
1. Еремин, М. А. Конечно-объемная схема интегрирования уравнений гидродинамики / М. А. Еремин, А. В. Хоперсков, С. А. Хоперсков // Изв. Волгогр. гос. техн. ун-та. —
2010. — № 6:8. — C. 24-27.
2. Жумалиев, А. Г. Численная схема cSPH-TVD: моделирование фронта ударной волны / А. Г. Жумалиев, С. С. Храпов // Вестник Волгоградского государственного университета. Серия 1, Математика. Физика. — 2012. — Т. 17, № 2. — C. 60-67.
3. Кузьмин, Н. М. Численная схема cSPH-TVD: исследование влияния ограничителей наклонов / Н. М. Кузьмин, А. В. Белоусов, Т. С. Шушкевич, С. С. Храпов // Вестник Волгоградского государственного университета. Серия 1, Математика. Физика. — 2014. — № 1. — C. 22-33.
4. Макарян, В. Г. Слабые ударные волны в неравновесных средах с отрицательной дисперсией / В. Г. Макарян, Н. Е. Молевич // Журнал технической физики. — 2005. — Т. 75, вып. 6. — C. 13-18.
5. Молевич, Н. Е. Волны в среде с отрицательной второй вязкостью / Н. Е. Молевич, А. Н. Ораевский // Труды ФИАН СССР. — 1992. — Т. 222. — C. 45-95.
6. Осипов, А. И. Неравновесный газ / А. И. Осипов // Соросовский образовательный журнал. — 1998. — № 7. — C. 95-101.
7. Храпов, C. C. Численная схема для моделирования динамики поверхностных вод на основе комбинированного SPH-TVD-подхода / C. C. Храпов, А. В. Хоперсков, Н. М. Кузьмин, А. В. Писарев, И. А. Кобелев // Вычислительные методы и программирование. —
2011. — Т. 12, № 1. — C. 282-297.
8. Шушкевич, К. С. Одномерная численная схема для газодинамического моделирования на основе комбинированного подхода SPH-PPM / К. С. Шушкевич, Н. М. Кузьмин // Вестник магистратуры. — 2013. — № 5 (20). — C. 40-44.
9. Leer, B. van. Towards the ultimative conservative difference scheme. III. Upstream-centered finite-difference schemes for ideal compressible flow / B. van. Leer // Journal of Computational Physics. — 1977. — Vol. 23, № 3. — P. 263-275.
10. Roe, P. L. Efficient construction and use of approximate Riemann solvers / P. L. Roe, J. Pike // Computing Methods in Applied Sciences and Engineering. — Amsterdam : North-Holland, 1983. — Vol. VI. — P. 499-518.
11. Toro, E. F. Riemann solvers and numerical methods for fluid dynamics / E. F. Toro. — Berlin ; Heidelberg : Springer-Verlag, 1999. — 624 p.
12. Toro, E. F. Godunov Methods: Theory and Applications / E. F. Toro. — New York : Kluwer Academic/Plenum Publ., 2012. — 1077 p.
REFERENCES
1. Eremin M.A., Khoperskov A.V., Khoperskov S.A. Konechno-obyemnaya skhema integrirovaniya uravneniy gidrodinamiki [Finite volume sheme of integration for hydrodynamics equations]. Izv. Volgogr. gos. tеkhn. un-ta, 2010, no. 6:8, pp. 24-27.
2. Zhumaliev A.G., Khrapov S.S. Chislennaya skhema cSPH-TVD: modelirovanie fronta udarnoy volny [Numerical scheme based on the combined SPH-TVD appoach: simulation of the shock front]. Vеstnik Volgogradskogo gosudars^nnogo univеrsitеta. Sеriya 1, Matеmatika. Fizika [Science Journal of Volgograd State University. Mathematics. Physics], 2012, vol. 17, no. 2, pp. 60-67.
3. Kuzmin N.M., Belousov A.V., Shushkevich T.S., Khrapov S.S. Chislennaya skhema cSPH-TVD: issledovanie vliyaniya ogranichiteley naklonov [Numerical scheme CSPH-TVD: investigation of influence slope limiters]. Vеstnik Volgogradskogo gosudars^nnogo univеrsitеta. Sеriya 1, Matеmatika. Fizika [Science Journal of Volgograd State University. Mathematics. Physics], 2014, no. 1, pp. 22-33.
4. Makaryan V.G., Molevich N.E. Slabye udarnye volny v neravnovesnykh sredakh s otritsatelnoy dispersiey [Weak shock waves in nonequilibrium media with negative dispersion]. Zhurnal tеkhnichеskoy fiziki [Technical Physics], 2005, vol. 75, iss. 6, pp. 13-18.
5. Molevich N.E., Oraevskiy A.N. Volny v srede s otritsatelnoy vtoroy vyazkostyu [Waves in a medium with a negative second viscosity]. Trudy FIAN SSSR [Proceedings of FIAN the USSR], 1992, vol. 222, pp. 45-95.
6. Osipov A.I. Neravnovesnyy gaz [Nonequilibrium gas]. Sorosovskiy obrazovatеlnyy zhurnal, 1998, no. 7, pp. 95-101.
7. Khrapov C.C., Khoperskov A.V., Kuzmin N.M., Pisarev A.V., Kobelev I.A. Chislennaya skhema dlya modelirovaniya dinamiki poverkhnostnykh vod na osnove kombinirovannogo SPH-TVD-podkhoda [Numerical scheme for the simulation of the dynamics of surface water based on the combined SPH-TVD-approach]. Vychislitеlnyе mеtody i programmirovaniе, 2011, vol. 12, no. 1, pp. 282-297.
8. Shushkevich K.S., Kuzmin N.M. Odnomernaya chislennaya skhema dlya gazodinamicheskogo modelirovaniya na osnove kombinirovannogo podkhoda SPH-PPM [One-dimensional numerical scheme for gas-dynamical simulations based on the combined approach SPH-PPM]. Vеstnik magistratury, 2013, no. 5 (20), pp. 40-44.
9. Leer B. van. Towards the ultimative conservative difference scheme. III. Upstream-centered finite-difference schemes for ideal compressible flow. Journal of Computational Physics, 1977, vol. 23, no. 3, pp. 263-275.
10. Roe P.L., Pike J. Efficient construction and use of approximate Riemann solvers. Computing Methods in Applied Sciences and Engineering, Amsterdam, North-Holland Publ., 1983, vol. VI, pp. 499-518.
11. Toro E.F. Riemann solvers and numerical methods for fluid dynamics. Berlin ; Heidelberg, Springer-Verlag, 1999. 624 p.
12. Toro E.F. Godunov Methods: Theory and Applications. New York, Kluwer Academic/Plenum Publ., 2012. 1077 p.
NUMERICAL MODELING OF THE DYNAMICS OF SOUND WAVES IN ACTIVE ENVIRONMENTS USING THE MUSCL SCHEME
Bochka^va Ekatеrina Viktorovna
Student, Institute of Mathematics and Information Technologies, Volgograd State University [email protected], [email protected] Prosp. Universitetsky, 100, 400062 Volgograd, Russian Federation
Khrapov Sеrgеy Sеrgееvich
Candidate of Physical and Mathematical Sciences, Associate Professor, Department of Information Systems and Computer Modeling, Volgograd State University [email protected], [email protected]
Prosp. Universitetsky, 100, 400062 Volgograd, Russian Federation
Abstract. The results of the Godunov type method using numerical schemes of the first and second level of precision are compared. Influence of slope limiters and flux computation methods to quality of numerical solution are investigated.
Four versions of slope limiters are investigated: minmod, van Leer, van Albada, superbee. Three methods of numerical flux computation also investigated: Lax — Friedrichs, Harten — Lax — van Leer and Harten — Lax — van Leer — Contact.
It is shown, that superbee leads to very similar numerical solution quality. Also, the results show, that the use of the MUSCL scheme gives smaller error calculations in 2 times, regardless of the method numerical solution of Riemann problem.
The instability of acoustic waves in active environments are investigated on the basis of the realisation of the method HLLC and MUSCL schemes.The results indicate significant effects of stationary nonequilibrium (generated through its new viscosity-dispersive and nonlinear properties of the medium) on the evolution of the gas-dynamic perturbations.
Key words: numerical schemes, slope limiters, method of Godunov type, MUSCL, HLLC.
