Особенности расчета взаимодействия ударной волны с газовым пузырем в мелкодисперсной газовзвеси

2019 Математика и механика № 57

УДК 533.6.011.72; 519.63 Б01 10.17223/19988621/57/8

Д.В. Садин, В.А. Давидчук

ОСОБЕННОСТИ РАСЧЕТА ВЗАИМОДЕЙСТВИЯ УДАРНОЙ ВОЛНЫ С ГАЗОВЫМ ПУЗЫРЕМ В МЕЛКОДИСПЕРСНОЙ ГАЗОВЗВЕСИ

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

Ключевые слова: ударная волна, газовзвесь, газовый пузырь, неустойчивость Рихтмайера - Мешкова, разностная схема

Потоки газовзвесей находят широкое применение в практических приложениях [1-3]. Для таких явлений характерно наличие существенного диапазона изменения газодинамических параметров течений, границы раздела сред, неоднород-ностей по плотности и крупных включений. Несомненный интерес представляют структурно-сложные задачи, которые отличают рефракция ударных волн на границе раздела сред и одновременно развитие неустойчивости и турбулентности.

При численном моделировании ударно-волновых процессов в газовзвесях с мелкодисперсной фракцией возникает ряд проблем. Для таких течений характерно явление жесткости, связанное с существенным различием времен релаксации фаз. Применение традиционных явных алгоритмов расчета становится практически не пригодным ввиду неоправданно малого шага по времени. Для преодоления этой трудности более 20 лет назад предложены неявные и полунеявные методы расчета источниковых членов (межфазных взаимодействий) [4-6]. На модельных односкоростных и двухскоростных уравнениях в частных производных показано, что такой подход справедлив для широкого класса разностных схем [4, 5]. Продемонстрирована высокая устойчивость схем даже для некорректно поставленных разностных задач - возможность устойчивого проведения расчетов с числом Куранта больше единицы СБЬ>1. В частности, для линейных задач построены аппроксимирующие схемы, доказана их устойчивость при СБЬ > 1, но из этих двух свойств не следует сходимость [6]. Проведенный анализ схем для двухскоростной динамики гетерогенных сред показал, что для практически важных случаев достаточно неявно рассчитывать только быстрые компоненты решения (газовую фазу). Изучены влияние на допустимый шаг по времени (число Куранта) шага расчетной сетки, размеров частиц газовзвеси, концентрации дисперсной фазы и других величин [5]. Введено понятие К-устойчивости схемы - условия устойчивости не зависят от интенсивности межфазных взаимодействий и размеров сетки [7, 8]. Другой подход основан на расщеплении законов сохранения газовзвесей на однородную систему уравнений (без источников) и введение корректирующего шага для

релаксационных процессов в виде обыкновенных дифференциальных уравнений [9]. Однако в работе [10] показано, что указанная стратегия построения разностной схемы может приводить к существенному дисбалансу градиентно-потоковых и источниковых членов уравнений, например, при получении численного решения методом установления. Интерес в исследованиях при численном моделировании ударно-волновых течений газовзвесей вызывает явление жесткости [11].

Весьма важным требованием, предъявляемым к разностным схемам, является возможность разрешения деталей ударной волновой картины и часто сопровождаемых явлений развития неустойчивости на контактных поверхностях и турбулентности. Схемы первого порядка точности являются избыточно диссипативны-ми и их применение требует значительных временных затрат для получения требуемого качества решения для структурно-сложных течений газа и газовзвесей. Для сохранения монотонности при повышении порядка аппроксимации схемы используют схемы высокой разрешающей способности типа Годунова [9, 12], ТУи [10], WENO-типа [13], с адаптивной искусственной вязкостью [14] и др.

Еще ряд проблем связан с численным моделированием взаимодействия волн с границами раздела сред с различными уравнениями состояния и сходимостью численных решений к равновесным течениям при уменьшении размеров частиц газовзвеси. Газодинамические течения в рамках сплошной среды в общем случае описываются уравнениями Навье-Стокса. Численное моделирование таких течений является корректным при условии — численная диффузия существенно меньше физической. В настоящей работе используется подход [13, 15] в рамках уравнения Эйлера, например для оценки разрешения сеток, обеспечивающих в последующих исследованиях корректное моделирование на основе уравнений Навье-Стокса при заданном числе Рейнольдса [13].

Постановка задачи

Для описания ударно-волновых процессов в газовзвесях запишем уравнения двухскоростной двухтемпературной динамики двухфазной среды [16] при известных допущениях в виде [14]

^е+в

Ч = [ Р2 , Р1^1, Р2^ Р2е2 , Р1е1 + Р2к2 ]т ; 6 = [ ^ Р2Г2 , Р1 ^ V, Р2Г2Г2 , Р2е2Г2 , Р1е1V + Р2к2Г2 ]т ^ Г = [0,0, р, р, 0, р (V + а2У2 ) , н = [0,0, — г, Г, дт, —ег ]т ;

= Шая (V- , V-, V, V, V., V-), в = Шая [1,1, а1, а 2,1,1];

р. = р°аг (/ = 1,2), Е1 = е1 + /2, К2 = у^/2 .

Здесь и далее индексы 1 и 2 внизу относятся соответственно к параметрам несущей и дисперсной фаз, индекс ° сверху относится к истинным значениям плотности; V - оператор Гамильтона. Через а., р., V., Ei, е., р обозначены объемная доля, приведенная плотность, вектор скорости, полная и внутренняя энергии единицы массы /-й фазы, давление газа; Г, - соответственно вязкая составляющая

силы межфазного взаимодействия, мощность теплообмена между газом и частицами в единице объема; / — время.

Для замыкания системы (1) используем уравнения состояния идеального кало-рически совершенного газа и несжимаемых твердых частиц: p = (у1 - 1)р1е1,

e1 = cv71, e2 = c2T2, (yj, cv, c2, p2} = const, где T1, T2 - температура несущей фазы и частиц; yj; cv - показатель адиабаты и удельная теплоемкость газа при постоянном объеме; c2 - удельная теплоемкость частиц. Интенсивности межфазного трения и теплообмена F^, QT задаются на основе экспериментальных соотношений

[17-19]. Начальные условия задаются применительно к рассматриваемым ниже задачам. Граничные условия: на стенках - отражения, на внешних границах (во всех рассматриваемых ниже задачах слева и справа расчетной области) - экстраполяция.

Метод решения

Для численного решения используем схему с настраиваемыми диссипативны-ми свойствами [14, 20] в конечно-объемной реализации [21] с аппроксимацией второго порядка по времени и пространству на гладких решениях — CDP2 (customizable dissipative properties).

Переход на следующий временной слой выполняется двухшаговым TVD-методом Рунге-Кутты [22]:

q(1) = qk + tL (qk ) , qk+1 = 0.5 (qk + q(1)) + 0.5tL (q(1)) .

При этом на каждом шаге производится расщепление пространственного разностного оператора L (q) по физическим процессам с использованием искусственной вязкости Q^1/2 типа Христенсена [23] с ограничителями вязкости yv [20, 21]:

1) qk+1/2 -SB(qkn+V2)т = qkn +(1 -S)B(qk„ )т-Бкп (Fkn+V2-Fkn_V2)т/ , где ^±1/2 =[0,0, p)kn±j/2, p n±1/ 2, 0, p n±1/ 2 KV +a 2v2 )П±1/2 ^ ,

pn±1/2 = pn±1/2 + Qn±1/2 ;

2) q<n') = qk+1/2-x[Vd G ]]+-.

Здесь т, h - временные и пространственные шаги сетки; верхний индекс k относится к временному слою; нижний целый индекс n относится к центру ячейки, а полуцелые n±1/2 к ее граням; [VdG]k+12 - TVD-реконструкция потоков путем

взвешенной линейной комбинации противопоточной и центральной аппроксимаций конвективных членов с ограничителями потоков у f [21]; S- весовой множитель источниковых членов (S = 1 неявный учет межфазных взаимодействий, S = 0 - явный).

Высокая устойчивость обеспечивается учетом межфазных взаимодействий на верхнем временном слое (S = 1), что ранее было предложено и обосновано в работах [4, 5, 14 и др.]. Искомые функции рассчитываются явно (без итераций) при линейной зависимости источниковых членов от решения. Однако при использо-

вании эмпирических критериальных соотношений такие зависимости - нелинейные. Например, сила межфазного трения определяется суммой слагаемых, зависящих от разности скоростей фаз в различных степенях. Для исключения итерационных процедур используем линеаризацию источников, учитывая неявно линейную часть [5]. Как показал опыт расчетов течений газовзвесей в широком диапазоне параметров потоков и размеров частиц (от 100 нм и выше), запас устойчивости алгоритма при таком подходе сохраняется. Число Куранта для всех задач в настоящей работе принято CFL = 0.4 , а ограничители искусственной вязкости и потоков - Van Leer [21].

Верификация метода

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

Первая из них - задача распада разрыва в газовзвеси, имеющая точное автомодельное решение в равновесном приближении [24, 14]. Тестовая задача решается на интервале x / Ц е (0,1) с начальным разрывом в точке x0 / Ц = 0.48 и на отрезке времени t е [0, t1 ], где Ц = 1 м - размер расчетной сетки с разрешением h / Ц = 1/500 ; tj = 1.2 мс - время окончания счета. Начальные условия задачи в системе СИ представлены в табл. 1. Газодинамические величины помечены индексом L - слева от разрыва, а индексом R - справа. Двухфазная смесь в начальный момент находится в состоянии покоя. Частицы - монодисперсные, диаметром d = 0.1 мкм , поровый газ - воздух с показателем адиабаты у = 7 / 5.

Таблица 1

Начальные условия тестовой задачи распада разрыва в газовзвеси

a2L PlL Pl a2R PlR Pr

0.001 11.8919 106 0.001 1.18919 105

После распада начального разрыва влево по газовзвеси распространяется волна разрежения г (рис. 1), а вправо ударная волна s, разделенные комбинированным разрывом (скачком пористости) с и зоной постоянного течения между ними. Для заданных мелких частиц времена скоростной и тепловой релаксаций малы, много меньше шага по времени, поэтому течение практически равновесное и сходится к точному решению (рис. 1, сплошные кривые).

Второй одномерный тест представляет собой модификацию известной задачи Карни и Кёрка [25]. В нашем случае ударная волна распространяется по газовзвеси и взаимодействует с ограниченной областью «чистого» газа (воздуха). Особенностью данной задачи является прохождение ударной волны и ее взаимодействие со средами, описываемыми различными уравнениями состояний. Кроме того, вычислительные трудности при решении этой задачи связаны с большим перепадом плотности.

x / l x / L

Рис. 1. Сравнение численных решений (пунктирные кривые) с точными (сплошные линии) задачи о распаде разрыва в газовзвеси в момент времени tf = 1.2 мс: a — величины давления 1 и плотности смеси 2, отнесенные к соответствующим начальным параметрам слева от разрыва; b — скорость смеси, отнесенная к скорости звука в газе справа от разрыва Fig. 1. Comparison of numerical (dashed lines) and exact solutions (solid lines) to the problem of decay of discontinuity in a gas suspension at the time instant of tf=1.2 ms: (a) 1, the pressure and 2, the density of mixture which are related to the corresponding initial parameters to the left of discontinuity; (b) the velocity of mixture related to a sound velocity in a gas to the right of discontinuity

В начальный момент времени ударная волна, удовлетворяющая соотношениям Рэнкина - Гюгонио в равновесной газовзвеси [14], расположена в точке x1 /L2 = 0.25 , а газовая область между x2 /L2 = 0.4 и x3 /L2 = 0.6 , где L2 = 1 м -

расчетная область. Как и в предыдущей задаче частицы имеют размер d = 0.1 мкм , а поровый газ — воздух. Исходные данные в системе СИ сведены в табл. 2.

Расчет выполнялся до момента времени t2 = 2 мс на сетке с разрешением h / L2 = 1/800 . Сравнение производилось с решением на достаточно мелкой сетке h / L2 = 1/4000, принятым за эталон. Результаты представлены на рис. 2 в виде распределений относительной плотности и скорости плотности смеси. Ко времени окончания счета реализуются пять распадов разрыва на контактных поверхностях с образованием волны разрежения и двух ударных волн слева от газовой области и двух ударных волн — справа. Полученное численное решение сохраняет монотонность с хорошим разрешением указанных деталей течения.

Таблица 2

Начальные условия модифицированной задачи Карни и Кёрка

Характерные области a2 Pi u P

x/L2 < 0.25 2.21214-10-3 2.63384 144.163 247467

0.25 < x/L2 < 0.4 и x/L2 > 0,6 10-3 1.18919 0 100000

0.4 < x/L2 < 0.6 10-10 1.18919 0 100000

1.0

0.8

Р / Р,

s

\

0.55

0.50

u / a„

0.6

0.45

0.4

0.40

0.2

0.4

0.6

0.8

1.0

0.35

0.2

0.4

0.6

0.8

1.0

* / ¿2 х / L2

Рис. 2. Сравнение численных решений (пунктирные кривые) с эталонными (сплошные линии) модифицированной задачи Карни и Кёрка в газовзвеси в момент времени t2 = 2 мс: a -плотность смеси, отнесенная к этой величине за ударной волной; b - скорость смеси, отнесенная к скорости звука в газе перед ударной волной

Fig. 2. Comparison of numerical (dashed lines) and benchmark solutions (solid lines) to a modified Kami and Quirk problem in the case of gas suspension at the time instant of t2 = 2 ms: (a) the density of mixture related to that after shock wave and (b) the velocity of mixture related to a sound velocity in a gas before shock wave

b

a

Еще одна важная особенность связана с пространственными эффектами интерференции волн разрежения, ударных волн и контактных разрывов, разделяющих два газа с различными показателями адиабаты. Рассмотрим показательную с этой точки зрения задачу о взаимодействии ударной волны с цилиндрическим пузырем гелия, которая нашла широкое применение при тестировании разностных схем. Ударная волна с числом Маха 1.22 движется влево по каналу с поперечным размером 8.9 см и встречает пузырь диаметром 5 см смеси гелия и воздуха с массовой долей 28 %. Перед ударной волной воздух и воздушно-гелиевая смесь находятся в состоянии покоя при температуре 293.23 К и давлении 101325 Па.

Результаты расчета по схеме СБР2 в виде численных шлирен-изображений функции градиента плотности газа показаны на рис. 3, а, Ь, с. Время отсчитывает-ся от момента прихода ударной волны на границу пузыря. Для сравнения приведем численные данные, полученные по схеме ЖБМО-типа четвертого порядка точности по пространству и третьего порядка по времени [15] (рис. 3, ё, е, /). Расчеты выполнены при одинаковых исходных данных, числе Куранта и дискретизации сеток с шагом Н=50 мкм. Численные решения хорошо согласуются между собой. Схема СБР2 демонстрирует возможности по разрешению тонких деталей ударно-волновой структуры, обладает малой диссипацией при сохранении монотонности решения. При этом вычислительные затраты существенно меньше, чем при использовании цитируемой ЖБМО-схемы.

Результаты расчета

Продолжим рассмотрение задачи распада разрыва в газовзвеси при изменении размеров частиц, когда эффекты неравновесности становятся существенными. Задача решается при исходных данных (см. табл. 1), а диаметры частиц заданы 10 и 100 мкм. После распада разрыва возникает релаксирующее течение газовзвеси.

Рис. 3. Численные шлирен-изображения функции градиента плотности газа в моменты времени: a, d - 62; b, e - 80; с, f - 245 мкс. Верхний ряд (a-c) - схема CDP2, нижний ряд (d-f) - WENO из [15]

Fig. 3. Numerical schlieren-images of the function of gas density gradient at the time instants of (a),(d) 62; (b),(e) 80; and (c),(f) 245 ^s. Upper row (a-c) corresponds to the CDP2 scheme; bottom row (d-f), the WENO from [15]

Вправо по поровому газу движется ударная волна уменьшающейся интенсивности (рис. 4, s). Влево распространяется волна разрежения (рис. 4, r). К моменту времени t2 =2 мс ударная волна в газовзвеси с частицами d = 100 мкм имеет скачок конечной интенсивности (рис. 4, штриховые кривые). Комбинированный разрыв с, который образовывался в равновесной газовзвеси (рис. 1, a), распадается на контактный разрыв в поровом газе i и скачок пористости с (рис. 4, a). В смеси с частицами d=10 мкм ударная волна вырождается в звуковую волну (рис. 4, штрихпунктирные кривые). При уменьшении размеров частиц зоны релаксации монотонно уменьшаются, а численное решение сходится к точному (рис. 4, сплошные линии).

В практически важных случаях ударно-волновые возмущения распространяются в газовзвесях, неоднородных по плотности. Поэтому задача взаимодействия ударной волны с газовым пузырем представляет как теоретический интерес, так и прикладное значение. Геометрические параметры соответствуют задаче взаимодействия ударной волны с гелиевым пузырем. А начальные условия - см. табл. 2 (за исключением скорости, которая имеет обратный знак). Задача решалась на равномерной сетке 4500x534 ячеек с одинаковыми шагами по пространственным координатам.

Плоская ударная волна распространяется в канале по мелкодисперсной (d = 0.1 мкм) газовзвеси и в некоторый момент времени встречает на своем пути неоднородность в виде воздушного пузыря. При этом возникает распад разрыва: с ударной волной внутри пузыря и ее рефракцией во внешнюю область, а также отраженной волной разрежения в противоположном направлении (см. рис. 5, а).

x / A x / L

Рис. 4. Сравнение численных решений

Рис. 5. Численные шлирен-изображения функции градиента плотности дисперсной фазы в моменты времени: a - 180; b - 250; с - 450; d - 650; e - 850; f - 1300 мкс Fig. 5. Numerical schlieren-images of the function of density gradient for a dispersed phase at the time instants of (a) 180; (b) 250; (с) 450; (d) 650; (e) 850; and f) 1300 jis

Скорость звука в газе выше, чем в газовзвеси. Поэтому ударная волна внутри пузыря движется быстрее, чем в окружающей газодисперсной среде. В момент времени от начала расчета /=250 мкс (рис. 5, b) взаимодействующая ударная волна прошла цилиндрическую неоднородность с образованием тройной точки. На рис. 5, c отмечается начало зарождающей неустойчивости Рихтмайера - Мешкова на поверхности комбинированного разрыва. Последующие моменты времени характеризуются более сложной ударно-волновой картиной и развитием крупномасштабной турбулентности, что хорошо видно на приведенных шлирен-изобра-жениях (рис. 5, d-f).

Выводы

Исследование проведено в рамках модели двухскоростной двухтемпературной динамики газовзвеси. Для расчета использована разностная схема со вторым порядком точности по пространству и времени CDP2, пригодная для решения жестких задач. Высокая устойчивость схемы обеспечена неявным учетом межфазных взаимодействий. Схема верифицирована на тестовой задаче распада разрыва в газовзвеси в сравнении с точным автомодельным решением в равновесной газовзвеси, модифицированной задаче Карни и Кёрка, а также сопоставлены с данными других авторов по взаимодействию ударной волны с цилиндрической неоднородностью в виде гелиевого пузыря. Тесты подтвердили малую диссипацию и монотонность численных решений, полученных по схеме CDP2. Рассмотрены релакси-рующие течения газовзвесей при изменении размеров частиц и их сходимость к равновесным решениям. Численно решена задача о распространении ударной волны в мелкодисперсной газовзвеси и ее взаимодействии с пузырем воздуха. Рассмотрены ударно-волновая картина, развитие неустойчивости Рихтмайера -Мешкова на поверхности раздела сред и образование крупномасштабной турбулентности. Проведенные исследования позволяют оценить численную диссипацию на заданных разрешениях сеток для корректного моделирования течений газовзвесей в рамках Навье - Стокса при заданном числе Рейнольдса, что является направлением дальнейших исследований.

ЛИТЕРАТУРА

1. Crowe C.T., Schwarzkopf J.D., Sommerfeld M., Tsuji Y. Multiphase flows with droplets and particles. New York: CRC Press, 2012. 487 p.

2. Алхимов А.П., Клинков С.В., Косарев В. Ф., Фомин В.М. Холодное газодинамическое напыление. Теория и практика / под ред. В.М. Фомина. М.: Физматлит, 2010. 536 с.

3. Gidaspow D. Multiphase Flow and Fluidization. New York: Academic Press, 1994. 467 p.

4. Садин Д.В. Модифицированный метод крупных частиц для расчета нестационарных течений газа в пористой среде // Ж. вычисл. матем. и матем. физ. 1996. Т. 36. № 10. С. 158-164.

5. Садин Д.В. Метод расчета волновых гетерогенных течений с интенсивным межфазным взаимодействием // Журн. вычисл. матем. и матем. физ. 1998. Т. 38. № 6. С. 1033-1039.

6. Садин Д.В. О сходимости одного класса разностных схем для уравнений нестационарного движения газа в дисперсной среде // Журн. вычисл. матем. и матем. физ. 1998. Т. 38. № 9. С. 1572-1577.

7. Садин Д.В. О жесткости систем дифференциальных уравнений в частных производных, описывающих движения гетерогенных сред // Математ. моделир. 2002. Т. 14. № 11. С. 43-53.

8. Садин Д.В. Проблема жесткости при моделировании волновых течений гетерогенных сред с трехтемпературной схемой межфазного тепло- и массообмена // Прикл. механ. и техн. физ. 2002. Т. 43. № 2. С. 136-141. DOI: 10.1023/A:1014714012032.

9. Saurel R. and Abgrall R. A multiphase Godunov method for compressible multifluid and multiphase flows // J. Comput. Phys. 1999. 150(2). P. 425-467. DOI: 10.1006/jcph. 1999.6187.

10. Gascon Ll., Corberan J.M. Construction of Second-Order TVD Schemes for Nonhomogene-ous Hyperbolic Conservation Laws // J. Comput. Phys. 2001. 172. P. 261-297. DOI: 10.1006/jcph.2001.6823.

11. Миньков Л.Л., Гольдина Н.В. Особенности численного решения задачи о распространении ударной волны по газовзвеси с мелкими частицами // Вестн Том. гос. ун-та. Математика и механика. 2017. № 49. С. 94-104. DOI: 10.17223/19988621/49/9.

12. Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics. 3rd ed. Berlin: Springer-Verlag, 2009. 724 p. DOI: 10.1007/b79761.

13. Shi J., Zhang Y.-T., Shu C.-W. Resolution of high order WENO schemes for complicated flow structures // J. Comput. Phys. 2003. V. 186. P. 690-696. DOI: 10.1016/S0021-9991(03) 00094-9.

14. Садин Д.В. TVD-схема для жестких задач волновой динамики гетерогенных сред негиперболического неконсервативного типа // Журн. вычисл. матем. и матем. физ. 2016. Т 56. № 12. С. 2098-2109. DOI: 10.7868/S0044466916120152.

15. Coralic V. and Colonius T. Finite-volume WENO scheme for viscous compressible multi-component flows // J. Comput. Phys. 2014. V. 274. P. 95-121. DOI: 10.1016/j.jcp. 2014.06.003.

16. Нигматулин Р.И. Динамика многофазных сред. Ч. 1, 2. М.: Наука, 1987.

17. Ergun S. Fluid flow through packed columns // Chem. Eng. Progress. 1952. V. 48. No. 2. P. 89-94.

18. ЧудновскийА.Ф. Теплообмен в дисперсных средах. М.: Гостехтеориздат, 1954.

19. Стернин Л.Е., Маслов Б.П., Шрайбер А.А., Подвысоцкий А.М. Двухфазные моно- и полидисперсные течения газа с частицами. М.: Машиностроение, 1980.

20. Садин Д. В. Схемы с настраиваемыми диссипативными свойствами для численного моделирования течений газа и газовзвесей // Матем. моделирование. 2017. Т. 29. № 12. С. 89-104.

21. Садин Д.В. Применение схемы с настраиваемыми диссипативными свойствами к расчету течений газа с развитием неустойчивости на контактной границе // Научно-технический вестник информационных технологий, механики и оптики. 2018. Т. 18. № 1. С. 153-157. DOI: 10.17586/2226-1494-2018-18-1-153-157.

22. Gottlieb S. and Shu C.-W. Total variation diminishing Runge-Kutta schemes // Mathematics of Computation. 1998. V. 67. No. 221. P. 73-85. DOI: 10.1090/S0025-5718-98-00913-2.

23. Christensen R.B. Godunov Methods on a Staggered Mesh - An Improved Artificial Viscosity. Technical Report UCRL-JC-105269. 1990. 11 p.

24. Иванов А.С., Козлов В.В., Садин Д.В. Нестационарное истечение двухфазной дисперсной среды из цилиндрического канала конечных размеров в атмосферу // Изв. РАН. Механика жидкости и газа. 1996. № 3. С. 60-66. DOI: 10.1007/BF02030221.

25. Abgrall R. How to Prevent Pressure Oscillations in Multicomponent Flow Calculations: A Quasi Conservative Approach // J. Comput. Phys. 1996. V. 125. P. 150-160. DOI: 10.1006/jcph.1996.0085.

Статья поступила 03.10.2018 г.

Sadin D.V., Davidchuk V.A. (2019) DISTINCTIVE FEATURES OF CALCULATING

INTERACTION BETWEEN SHOCK WAVE AND GAS BUBBLE IN A FINELY DISPERSED

GAS SUSPENSION. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika

[Tomsk State University Journal of Mathematics and Mechanics]. 57. pp. 100-111

DOI 10.17223/19988621/57/8

Keywords: shock wave, gas suspension, gas bubble, Richtmyer-Meshkov instability, difference

scheme.

The problem of interaction between shock wave and gas bubble in a finely dispersed gas suspension is studied using a two-velocity two-temperature formulation. The numerical method, which is applicable for a simulation of multiphase flows governed by the stiff Euler equations, is utilized. Implementation of the scheme is split into two phases. The first uses the central differences of both deformation and gradient terms with Christensen-type artificial viscosity. The total variation diminishing (TVD)-type reconstructions of the fluxes are used in the second phase applying a weighted linear combination of upwind and central approximations of convective terms with flux limiters. The second-order TVD Runge-Kutta (RK) algorithm is employed to march the solution in time. A high stability is ensured by either implicit or semi-implicit calculating method for the source terms in the equations, which have been proposed and developed over the last decades. The properties of elaborated numerical method are verified by considering several challenging one- and two-dimensional test problems as compared to the exact self-similar equilibrium solutions and to the results of other authors. A convergence to the equilibrium solutions is confirmed at various particle sizes. The shock-wave pattern, the Richtmyer-Meshkov instability developing along the bubble interface, and the large-scale turbulence generation are studied.

SADIN Dmitriy Viktorovich (Doctor of Technical Science, Professor, Mozhaysky Military Space Academy, Saint-Petersburg, Russian Federation). E-mail: sadin@yandex.ru

DAVIDCHUK Viktor Aleksanlrovich (Mozhaysky Military Space Academy, Saint-Petersburg, Russian Federation). E-mail:david_lxii@mail.ru

REFERENCES

1. Crowe C.T., Schwarzkopf J.D., Sommerfeld M., Tsuji Y. (2012) Multiphase Flows with Droplets and Particles. New York: CRC Press.

2. Alkhimov A.P., Klinkov S.V., Kosarev V.F., Fomin V.M. (2010) Kholodnoe gazodi-namicheskoe napylenie. Teoriya ipraktika [Cold gas-dynamic spraying. Theory and practice]. Moscow: Fizmatlit.

3. Gidaspow D. (1994) Multiphase Flow and Fluidization. New York: Academic Press.

4. Sadin D.V. (1996) A modified large-particle method for calculating unsteady gas flows in a porous medium. Comp. Maths Math. Phys. 36(10). pp. 1453-1458.

5. Sadin D.V. (1998) A method for computing heterogeneous wave flows with intense phase interaction. Comp. Maths Math. Phys. 38(6). pp. 987-993.

6. Sadin D.V. (1998) On the convergence of a certain class of difference schemes for the equations of unsteady gas motion in a disperse medium. Comp. Maths Math. Phys. 38(9). pp. 1508-1513.

7. Sadin D.V. (2002) O zhestkosti sistem uravneniy v chastnykh proizvodnykh, opisyva-yushchikh dvizheniya geterogennykh sred [On stiff systems of partial differential equations for motion of heterogeneous media]. Matematicheskoe modelirovanie - Mathematical Models and Computer Simulations. 14(11). pp. 43-53.

8. Sadin D.V. (2002) Stiffness problem in modeling wave flows of heterogeneous media with a three-temperature scheme of interphase heat and mass transfer. Journal of Applied Mechanics and Technical Physics. 43(2). pp. 286-290. DOI: 10.1023/A:1014714012032.

9. Saurel R., Abgrall R. (1999) A multiphase Godunov method for compressible multifluid and multiphase flows. J. Comput. Phys. 150(2). pp. 425-467. DOI: 10.1006/jcph.1999.6187.

10. Gascon Ll., Corberan J.M. (2001) Construction of second-order TVD schemes for nonhomo-geneous hyperbolic conservation laws. J. Comput. Phys. 172(1). pp. 261-297. DOI: 10.1006/jcph.2001.6823.

11. Min'kov L.L., Gol'dina N.V. (2017) Osobennosti chislennogo resheniya zadachi o rasprostranenii udarnoy volny po gazovzvesi s melkimi chastitsami [Peculiarities of a numerical solution of the problem of shock wave propagation over a gas suspension with small particles]. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika - Tomsk

State University Journal of Mathematics and Mechanics. 49. pp. 94-104. DOI: 10.17223/ 19988621/49/9.

12. Toro E.F. (2009) Riemann Solvers and Numerical Methods for Fluid Dynamics. 3rd ed. Berlin: Springer-Verlag. DOI: 10.1007/b79761.

13. Shi J., Zhang Y.T., Shu C.W. (2003) Resolution of high order WENO schemes for complicated flow structures. J. Comput. Phys. 186. pp. 690-696. DOI: 10.1016/S0021-9991(03) 00094-9.

14. Sadin D.V. (2016) TVD scheme for stiff problems of wave dynamics of heterogeneous media of nonhyperbolic nonconservative type. Comp. Maths Math. Phys. 56(12). pp. 2068-2078. DOI: 10.7868/S0044466916120152.

15. Coralic V., Colonius T. (2014) Finite-volume WENO scheme for viscous compressible multi-component flows. J. Comput. Phys. 274. pp. 95-121. DOI: 10.1016/j.jcp.2014.06.003.

16. Nigmatulin R.I. (1987) Dinamika mnogofaznykh sred [Dynamics of multiphase media]. Part. 1, 2. Moscow: Nauka.

17. Ergun S. (1952) Fluid flow through packed columns. Chem. Eng. Progress. 48(2). pp. 89-94.

18. Chudnovsky A.F. (1954) Teploobmen v dispersnykh sredakh [Heat transfer in dispersed media]. Moscow: Gostekhteorizdat.

19. Sternin L.E., Maslov B.P., Shrayber A.A., Podvysotskiy A.M. (1980) Dvukhfaznye mono- i polidispersnye techeniya gaza s chastitsami [Two-phase mono- and polydispersed flows of gas with particles]. Moscow: Mashinostroenie.

20. Sadin D.V. (2017) Skhemy s nastraivaemymi dissipativnymi svoystvami dlya chislennogo modelirovaniya techeniy gaza i gazovzvesey [Schemes with customizable dissipative properties as applied to gas-suspensions flow simulation]. Matematicheskoe modelirovanie -Mathematical Models and Computer Simulations. 29(12). pp. 89-104.

21. Sadin D.V. (2018) Application of scheme with customizable dissipative properties for gas flow calculation with interface instability evolution. Scientific and Technical Journal of Information Technologies, Mechanics and Optics. 18(1). pp. 153-157. DOI: 10.17586/22261494-2018-18-1-153-157.

22. Gottlieb S., Shu C.W. (1998) Total variation diminishing Runge-Kutta schemes. Mathematics of Computation. 67(221). pp. 73-85. DOI: 10.1090/S0025-5718-98-00913-2.

23. Christensen R.B. (1990) Godunov Methods on a Staggered Mesh -An Improved Artificial Viscosity. Technical Report UCRL-JC-105269.

24. Ivanov A.S., Kozlov V.V., Sadin D.V. (1996) Unsteady flow of a two-phase disperse medium from a cylindrical channel of finite dimensions into the atmosphere. Fluid Dynamics. 31(3). pp. 386-391. DOI: 10.1007/BF02030221.

25. Abgrall R. (1996) How to prevent pressure oscillations in multicomponent flow calculations: a quasiconservative approach. J. Comput. Phys. 125. pp. 150-160. DOI: 10.1006/ jcph.1996.0085.

Received: October 3, 2018