CC BY
28
УДК 519.63+532.529.5
DOI: 10.14529/ mmp190209
МОДИФИКАЦИЯ МЕТОДА КРУПНЫХ ЧАСТИЦ ДО СХЕМЫ ВТОРОГО ПОРЯДКА ТОЧНОСТИ ПО ПРОСТРАНСТВУ И ВРЕМЕНИ ДЛЯ УДАРНО-ВОЛНОВЫХ ТЕЧЕНИЙ ГАЗОВЗВЕСИ
Д.В. Садин, Военно-космическая академия имени А.Ф. Можайского, г. Санкт-Петербург, Российская Федерация
Развивается предложенный нами ранее подход построения разностных схем для решения жестких задач ударно-волновых течений гетерогенных сред с использованием неявного безытерационного алгоритма расчета межфазных взаимодействий. Метод крупных частиц модифицирован до схемы второго порядка точности по времени и пространству на гладких решениях. На первом этапе используются центральные разности с искусственной вязкостью ТУБ типа. На втором - ТУБ-реконструкция путем взвешенной линейной комбинации противопоточной и центральной аппроксимаций с ограничителями потоков. Схема дополнена двухшаговым методом Рунге - Кутта по времени. Схема является К-устойчивой - шаг по времени не зависит от интенсивности межфазных взаимодействий, а определяется числом Куранта для однородной системы уравнений (без источниковых членов). На тестовых задачах подтверждена монотонность, малая диссипативность, высокая устойчивость схемы и сходимость численных результатов к точным автомодельным равновесным решениям в газовзвеси. Показаны возможности схемы для численного моделирования физической неустойчивости и турбулентности. Метод может быть рекомендован для структурно-сложных течений газовзвесей.
Ключевые слова: течение газовзвеси; жесткая задача; 'разностная схема; устойчивость; точность.
Введение
Течения газовзвесей находят широкое применение в пневмотранспорте, аддитивных технологиях нанесения покрытий, системах порошкового пожаротушения, а также при анализе ударно-волновых процессов в аварийных ситуациях на угольных шахтах, мукомольных производствах и др. Математическое моделирование указанных процессов является эффективным инструментом оценки достижимого уровня технических параметров технологий, а в случае чрезвычайных ситуаций - потенциальных масштабов их последствий.
Метод крупных частиц для расчета газодинамических течений был предложен О.М. Белоцерковским и Ю.М. Давыдовым в начале 70-х годов прошлого века [1] и быстро завоевал популярность среди исследователей. Метод отличается простотой, надежностью и физичностью дискретного описания исследуемых процессов. Для расчета нестационарных волновых процессов в многофазных средах была предложена модификация метода крупных частиц [2]. Модифицированный алгоритм нашел широкое применение для решения разлета порошкообразных сред, взаимодействия ударных волн со слоями газовзвеси, обтекания смеси газа с частицами твердых тел и в других случаях [3]. Вместе с тем, при расчете двухфазных течений с интенсивным межфазным взаимодействием, например с малым (микронным и субмикронным) размером дисперсных частиц более двадцати лет назад возникла проблема, связанная с необходимостью неоправданного занижения шага по времени из условий устойчивости алгоритма. Это явление объяснялось жесткостью указанного класса задач, в
которых присутствуют решения с существенно различными характерными временами (временами релаксации фаз), а решение разделяется на быстрые и медленные компоненты. Для обеспечения приемлемой устойчивости предложены методы, основанные на неявном или полунеявном учете межфазных взаимодействий [4-8]. Следует отметить работы, в которых затрагиваются вопросы построения разностных схем для решения жестких задач динамики гетерогенных сред [9-12].
При отмеченных преимуществах метод крупных частиц обладает лишь первым порядком аппроксимации по пространству и времени. В настоящее время практика выдвигает задачи, требующие детального разрешения структурно-сложных пространственных гетерогенных течением с возможностью выявления и анализа развития неустойчивости. В этих случаях эффективность (вычислительные затраты для обеспечения заданной точности) применения метода крупных частиц снижается.
Настоящая работа посвящена модификации метода крупных частиц с повышением порядка аппроксимации до второго по времени и пространству на гладких решениях, пригодной для решения жестких задач газовзвесей.
1. Математическая модель газовзвеси
Рассмотрим динамику газовзвеси в рамках взаимопроникающих континуумов двухскоростной двухтемпературной двухфазной среды [3] при известных допущениях в виде [7]:
^ + VdG + B(VdF) = H(q), (1)
q = [Pl,P2,PlVl,P2V2,P2e2,PlEl + P2K2]T ,
G = [PlVi,P2V2,PIVIVI,P2V2V2, Р2в2V2,P?E?VI + P2^2V2]T ,
F = [0, 0,p,p, 0,p (alVl + «2V2)]T , H = [0, 0, —FM, FM, Qt, —Qt]T ,
Vd = diag (V-, V-, V, V, V-, V-), B = diag [1,1, a?, «2,1,1],
Pi = P>i(i = 1, 2), a? + a = 1, El = el + v?/2, K = v2/2.
Здесь и далее индексы 1 и 2 внизу относятся соответственно к параметрам несущей и дисперсной фаз, индекс о сверху относится к истинным значениям плотности; V -оператор Гамильтона. Через ai,Pi, vi, Ei,ei,p обозначены объемная доля, приведенная плотность, вектор скорости, полная и внутренняя энергии единицы массы i-ой фазы, давление газа; FM, Qt -соответственно вязкая составляющая силы межфазного взаимодействия, мощность теплообмена между газом и частицами в единице объема; t - время.
Для замыкания системы (1) используем уравнения состояния идеального кало-рически совершенного газа и несжимаемых твердых частиц: p = (yi — 1)р£el, el = cvT?, e2 = c2T2, {y?, cv, c2, p£} = const, где T?, T2 - температура несущей фазы и частиц; YbCv - показатель адиабаты и удельная теплоемкость газа при постоянном объеме; c2 - удельная теплоемкость частиц. Интенсивности межфазного трения и теплообмена FM, Qt задаются на основе экспериментальных соотношений [3]. Начальные условия задаются применительно к рассматриваемым ниже задачам. Граничные условия: на стенках - отражения, на внешних границах (во всех рассматриваемых ниже задачах слева и справа расчетной области) - экстраполяция.
2. Модификация метода крупных частиц
Запишем модифицированный алгоритм метода крупных частиц (без ограничения общности в одномерном случае), помечая разностное решение целыми индексами п в центре ячеек, полуцелыми п ± 1/2 на ее гранях, а верхним индексом к - временной слой.
1. Эйлеров этап
яП+1/2 - ¿и (яП+1/2) т = яП + (1 - ¿) н (яП) г - вп (ЁП+1/2 - Ё^) т/Л, (2)
F n± 1/2
О, 0,^1/2^1/2, 0, pn± 1/2 (a1v1 + a2v2)n±1/2 , pn± 1/2 = pn±1/2 + Qn±1/2
T ' ' k
где t, h - шаги сетки по времени и пространству, Q±1/2 - искусственная вязкость.
Модификация на первом этапе заключается в использовании искусственной вязкости TVD-типа с ограничителем вязкости ф [8, 13]:
Qn+1/2 = - (1 - Ф'с(гп+1/2)) Bv^llPn+l/2pÍ,n+l/2 ' Kn+1 ~ > (3)
r . (vk,n - <п-0 / K,n+1 - <п) , если (vk;n+1 - (рП+1 - рП) > 0,
r n+1/2 = S ( k k k k А
l,n+2 - v k,n+1) / (vk,n+1 - <n) , иначе,
где - коэффициент искусственной вязкости (для всех решаемых ниже задач принят = 1).
Приведем ограничители, используемые ниже:
0, метод крупных частиц,
Ф (r) = 4 (r + |r|) / (1 + r) , Van Leer, (4)
max [min (2r, 1), min (r, 2), 0], Superbee.
На гладких решениях ф (r) ^ 1, следовательно искусственная вязкость Q±1/2 ^ 0, и численное решение на первом этапе будет иметь второй порядок по пространству O (h2).
Учет межфазных взаимодействий осуществляется на верхнем временном слое при 5 =1 или нижнем при 5 = 0 [4, 5]. При линейной зависимости источниковых членов H (q) от искомых функций, например от разности скоростей фаз, решение (2) при 5 =1 выписывается явно. Для эмпирических замыкающих соотношений [3] такие зависимости имеют степенной вид. Для исключения итерационных алгоритмов будем проводить линеаризацию и учитывать на верхнем временном слое линейную часть [5].
2. Доработка метода крупных частиц на лагранжевом и заключительном этапах заключается в TVD-реконструкции искомых переменных ^ = {p¿, v, E ,e2} путем взвешенной линейной комбинации противопоточной и центральной аппроксимаций с ограничителем потоков ф/ [13]:
~k+1/2 = ( 1/2, если V>+V2 > 0,
f i,n+1 /2 =1 -
ИHаЧе,
^n+1/2 = - ф/ (riIn+1/^) 1/2 + ф/ (r+!n+1/^ ^k>++/12/2
ri,n+1/2
r+ i,n+1/2
fc+1/2 Pi.n
fc+1/2 Pi,n-1
Wk+1/2 ri,n+2
Wk+1/2 ri,n+1
rk+1/2 rt,n+1
fc+1/2 ; Pi.n
' i,n+1/2
Wk+1/2 г i,n+1
fc+1/2 ' Pi.n
Реконструированные переменные на гранях ячеек помечены тильдой , а
параметры с полуцелым пространственным индексом вычисляют как среднее их значений в центрах контрольных объемов = + 1/2 ^ /2. Расчет искомых переменных осуществляется следующим образом:
m = fc+V2 + ( Mk+V2 M
pi,n = pi,n + ^Mi,n-1/2 — Mi
fc+1/2 i,n+1/2
/h Mk+1/2 = pk+1/2 vk+1/2 /h, Mi,n±1/2 = Pi,n±1/2vi,n±1/2T,
,(1)
Pfc+1/2 k+1/2 I /vk+1/2 Mk+1/2 vk+1/2 M'
Pi,n Vi,n + (Дп-1/2Mi,n-1/2 — n /2M
fc+1/2 д «-fc+1/2
fc+1/2 i,n+1/2 Mi,n+1/2
e(1) = e2,n =
/h /pk
fc+1/2 fc+1/2 p2,n e2,n
+ I Pk+1/2 Mk+1/2 — Pk+1/2 Mk+1/2 ) /h /Pfc+1 + I e2,n-1/2M2,n-1/2 e2.n+1/2M2.n+1/2 i /h /p2
~"2,n+1/2 2,n+1/2
y2,n
(5)
E(1) = e(l) + Л,(?Л2 /2 Д P(l) E(1) = P(l) E(1)
E2,n = e2,n + I v2,n ) /2, Др2,пE2,n = p2,nE2,n
fc+1/2 Efc+1/2
p2,n E2
2,n
E
(1) 1,n
Pl
fc+1/2 p fc+1/2
pk+?/2 д ' Ei,n — ДР:
(i) E(?) 2,nE2,n
+ £ (E
i=1
Sfc+1/2 M fc+1/2 — E fc+1/2 M fc+1/2 vi,n-1/2Mi,n-1/2 Ei,n+1/2 Mi,n+1/2
/h
/ P'
fc+1 1,n '
В целом после лагранжевого, эйлерова и заключительного этапов при использовании известных ограничителей TVD-типа, например Van Leer или Superbee (4), имеем суммарный первый порядок по времени и второй по пространству на гладких решениях O (т + h2). Если задать — = 0 , то алгоритм <возвращается> к классическому методу крупных частиц и его модификации с неявным учетом межфазных взаимодействий [5] с первым порядком аппроксимации O (т + h).
Для повышения порядка аппроксимации по времени применим двухшаговой TVD-метод Рунге - Кутта [14]:
q(1) = qk + tL (qk) , (6)
qk+? = 0, 5 (qk + q(1)) + 0, 5tL (q(1)) . (7)
Здесь первый шаг (6) с разностным оператором L (qfc) представляет собой сокращенную запись алгоритма (2) - (5), а второй шаг (7) обеспечивает суммарную аппроксимацию O (т2 + h2).
За описанной дискретной моделью закрепилось название - схема с настраиваемыми диссипативными свойствами CDP2 (customizable dissipative properties) [8, 13], которое и используется ниже в этой работе. Настройка уровня численной диссипации может осуществляться выбором ограничителей —v, —/ и коэффициента искусственной вязкости, другой смысл этого названия заключается в автоматической подстройке (адаптации) схемы к локальным свойствам решения. Например, при появлении локальным максимумов (минимумов) в численном решении ограничитель вязкости обращается в ноль —v (r) = 0 , а величина искусственной вязкости (3) достигает максимального значения.
Схема является K-устойчивой. Устойчивость определяются условиями Куранта -Фридрихса - Леви для однородной системы уравнений (без источниковых членов), а
2
n
шаг по времени не зависит от интенсивности межфазных взаимодействий и размеров сетки [6, 8]. Схема прошла детальное тестирование на задачах газовой динамики, например [8, 13, 15].
3. Валидация схемы
Для проверки и подтверждения свойств схемы СВР2 (точности, устойчивости, уровня диссипативности и др.) рассмотрим тестовые задачи. Одна из них - комплексная одномерная задача, которая имеет точное решение в равновесном приближении [16, 17]. В начальный момент времени среда покоится и занимает расчетную область Ь = 0, 8 м , которая разделена на три части: 0 < х < х1 - камера высокого давления с газовзвесью; х1 < х < х2 - «чистый> газ (воздух); х2 < х < Ь - слой газовзвеси. Частицы газовзвеси являются сферическими, монодисперсными диаметром d = 0,1 мкм. Расчет выполнялся до момента времени £ =1 мс на сетке Л/Ь = 1/200. Начальные условия в системе СИ представлены в табл. 1.
Таблица 1
Начальные условия комплексной одномерной задачи
Характерные области cu2 PÍ щ,и2 p
0 < x/L < 0,4 10"3 11,8919 0 106
0,4 < x/L < 0,75 10-io 1,18919 0 105
0,75 <x/L< 1 10"3 1,18919 0 105
После распада разрыва в точке x1 влево по газовзвеси распространяется волна разрежения r (рис. 1). Вправо по газу движется ударная волна, которая через некоторое время встречает слой газовзвеси c2 в точке x2. Возникает второй распад разрыва с образованием двух ударных волн. Одна из которых распространяется вглубь слоя смеси газовзвеси s1, а отраженная - в противоположном направлении. После взаимодействия с границей раздела сред c1 происходит еще один распад разрыва с двумя ударными волнами s2 и Окончательная волновая конфигурация показана на рис. 1 в виде распределений относительных величин плотности смеси р/р0 (а) и скорости смеси м/а1 (б). Здесь начальные значения р0 - плотности смеси в камере высокого давления, а1 - скорости звука в воздухе перед разрывом. Сплошная кривая - точное автомодельное решение и результаты расчетов: штрихпунктирная линия - методом крупных частиц, пунктирная и штриховые кривые - схема CDP2 с ограничителем потоков Superbee и Van Leer соответственно.
Сравнение численных результатов (рис. 1) показывает, что метод крупных частиц является существенно диссипативным. Особенно это заметно при численном размазывании контактных и комбинированных разрывов c1 и c2. Для данной сетки схема CDP2 дает хорошее согласие с точным решением с меньшей диссипацией при использовании ограничителя потоков Superbee и несколько большее сглаживание разрывов (Van Leer).
Сходимость численных результатов в норме L1 к точному решению для метода крупных частиц и схемы CDP2 с разными ограничителями на примере комплексной одномерной задачи показана в табл. 2. Как видно из таблицы, ошибки расчета по схеме CDP2 в указанной норме на сетке 1/200 меньше, чем численные ошибки методом крупных частиц на более мелкой сетке с разрешением 1/800.
1,0 0,8 0,6 0,4 0,2
Р / Р0
и / а.
4ci
0 0,2 0,4 0,6 0,8 1,0
х / L1
а)
1,0 0,8 0,6 0,4 0,2
/ f S3
/ S1
/ г
/
х / L.
б)
Рис. 1. Сравнение численных результатов с точным решением (сплошная кривая): штрихпунктирная линия - метод крупных частиц; пунктир - CDP2 с ограничителем Superbee; штриховая линия - CDP2 (Van Leer)
Таблица 2
Ошибки численных решений в норме Li
г
Метод решения сетка 1/200 сетка 1/400 сетка 1/800
метод крупных частиц CDP2-Van Leer CDP2-Superbee 0,0276 0,0128 0,0106 0,0186 0,0063 0,0044 0,0133 0,0040 0,0028
Вторая тестовая задача представляет собой распространение ударной волны по газовой среде в канале с поперечным размером H = 0, 2 м, в котором имеется цилиндрическая область газовзвеси диаметром D = 0,1 м. Ввиду симметрии задачи расчеты выполнялись в верхней полуплоскости в области размером 5D х D с началом системы координат на оси справа. Начальное положение ударной волны - — 0, 6D, а центра цилиндрической области газовзвеси — 2D. Расчеты выполнялись на сетке h/D = 1/400. Начальные условия заданы в табл. 3.
Таблица 3
Начальные условия задачи взаимодействия ударной волны с газовзвесью
Характерные области «г PÍ Ml, U2 V
за ударной волной перед ударной волной областв газовзвеси Ю"10 Ю"10 Ю-3 2,43115 1,18919 1,18919 -281,769 0 0 284816 100000 100000
Эта задача является аналогом известного газодинамического теста взаимодействия ударной волны в воздухе с пузырем газа И.22 [18]. Экспериментальные данные и численные результаты широко используются для проверки работоспособности разностных схем, в частности для анализа их диссипативных свойств [19]. Для этого явления характерно различие в скоростях звука для воздуха и рефрижераторного газа И.22. В нашем случае также скорость распространения малых возмущений в равновесной двухфазной среде меньше скорости звука в воздухе.
Результаты расчетов для различных моментов времени показаны на рис. 2 в виде шлирен-изображений градиента плотности смеси, с использованием техники [20].
На этом рисунке сверху от оси симметрии приведены численные данные методом крупных частиц, а снизу CDP2 с ограничителем потоков Van Leer.
После начала взаимодействия возникает следующая конфигурация (рис. 2, а). Падающая ударная волна распадается на отраженную s1 от границы раздела сред, прошедшую s2 внутрь цилиндрической области газовзвеси и огибающую волну s3. В следующий момент времени (рис. 2, б) огибающая волна s3, имеющая большую скорость распространения, прошла зону облака газовзвеси. Скачок s2 движется медленнее и еще расположен внутри двухфазной среды. Ударная волна s1 отражается от стенок канала s4. В более поздние моменты времени волновая структура становится более сложной с образованием вторичных скачков s5, (рис. 2, гид), интерференции и развитием неустойчивости на поверхности разрыва сред us (рис. 2, д и е). Ударно-волновая картина и вихреобразование в газодисперсной среде качественно согласуется с аналогичным экспериментально подтвержденным явлением в двухком-понентном газе воздух - R22 [18, 19]. Схема CDP2 обладает малым уровнем диссипа-тивных свойств и хорошо разрешает детали течения (рис. 2, снизу от оси симметрии). Метод крупных частиц существенно размазывает комбинированный разрыв c (скачок пористости) и не позволяет выявить развитие неустойчивости на его поверхности us (рис. 2, сверху от оси симметрии).
г) д) е)
Рис. 2. Шлирен-изображения взаимодействия ударной волны с цилиндрической областью газовзвеси(^ = 0,1 мкм) в моменты времени: а) 280; б) - 400; в) 450; г) 500; д) 550; е) 1000 мкс. Сверху от оси симметрии - расчет методом крупных частиц; снизу - СБР2
Влияние фактора неравновесности при различии скоростей и температур фаз в поставленной задаче демонстрируется на рис. 3, где приведены результаты расчетов взаимодействия ударной волны с цилиндрической областью газовзвеси с диаметром частиц d = 2 мкм (верхний ряд - а, б, в) и d =10 мкм (нижний ряд - г, д, е). С увеличением размера частиц дисперсная фаза с меньшей интенсивностью вовлекается несущим газом вследствие увеличения времени скоростной и тепловой релаксации. Это приводит к уменьшению вихреобразования в дисперсной фазе - частицы <не замечают> турбулентности малого масштаба в газе.
г) д) е)
Рис. 3. Шлирен-изображения взаимодействия ударной волны с цилиндрической областью газовзвеси (d = 2 мкм - верхний ряд, d =10 мкм - нижний ряд) в моменты времени: а), г) - 400; б), д) 550; в), е) - 1000 мкс
Заключение
В рамках модели неравновесной динамики газовзвеси рассмотрена модификация метода крупных частиц до схемы второго порядка точности по пространству и времени. Схему отличает высокая устойчивость алгоритма за счет неявного расчета межфазных взаимодействий, малый уровень численной диссипации, монотонность и возможность детального разрешения особенностей структурно-сложных течений газовзвесей. Валидация алгоритма проведена на одномерной комплексной задаче, имеющей точное автомодельное решение для равновесной газовзвеси. Выполнен анализ сходимости численных решений в сравнении с методом крупных частиц. Результаты решения двумерной тестовой задачи взаимодействия ударной волны с цилиндрической областью газовзвеси продемонстрировали возможности схемы для моделирования явлений развития неустойчивости на границе раздела сред.
Литература
1. Белоцерковский, О.М. Нестационарный метод <крупных частиц> для газодинамических расчетов / О.М. Белоцерковский, Ю.М. Давыдов // Журнал вычислительной математики и математической физики. - 1971. - Т. 11, № 1. - С. 182-207.
2. Губайдуллин, А.А. Модифицированный метод крупных частиц для расчета нестационарных волновых процессов в многофазных средах / А.А. Губайдуллин, А.И. Ивандаев, Р.И. Нигматулин // Журнал вычислительной математики и математической физики. - 1977. - Т. 17, № 6. - С. 1531-1544.
3. Нигматулин, Р.И. Динамика многофазных сред / Р.И. Нигматулин. - М.: Наука, 1987.
4. Садин, Д.В. Модифицированный метод крупных частиц для расчета нестационарных течений газа в пористой среде / Д.В. Садин // Журнал вычислительной математики и математической физики. - 1996. - Т. 36, № 10. - С. 158-164.
5. Садин, Д.В. Метод расчета волновых гетерогенных течений с интенсивным межфазным взаимодействием / Д.В. Садин // Журнал вычислительной математики и математической физики. - 1998. - Т. 38, № 6. - С. 1033-1039.
6. Садин, Д.В. Проблема жесткости при моделировании волновых течений гетерогенных сред с трехтемпературной схемой межфазного тепло- и массообмена / Д.В. Садин // Прикладная механика и техническая физика. - 2002. - Т. 43, № 2. - С. 136-141.
7. Садин, Д.В. TVD-схема для жестких задач волновой динамики гетерогенных сред негиперболического неконсервативного типа / Д.В. Садин // Журнал вычислительной математики и математической физики. - 2016. - Т. 56, № 12. - С. 2098-2109.
8. Садин, Д.В. Схемы с настраиваемыми диссипативными свойствами для численного моделирования течений газа и газовзвесей / Д.В. Садин // Математическое моделирование. - 2017. - Т. 29, № 12. - С. 89-104.
9. Saurel, R. A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows / R. Saurel, R. Abgrall // Journal of Computational Physics. - 1999. - V. 150, № 2. -P. 425-467.
10. Saurel, R. Shock Jump Relations for Multiphase Mixtures with Stiff Mechanical Relaxation / R. Saurel, O. Le Metayer // Shock Waves. - 2007. - V. 16, № 3. - P. 209-232.
11. Saurel, R. Simple and Efficient Relaxation Methods for Interfaces Separating Compressible Fluids, Cavitating Flows and Shocks in Multiphase Mixtures / R. Saurel, F. Petitpas, R.A. Berry // Journal of Computational Physics. - 2009. - V. 228, № 5. - P. 1678-1712.
12. Грищенко, Д.С. Модификация метода крупных частиц для исследований течений газовзвесей / Д.С. Грищенко, Ю.М. Ковалев, Е.А. Ковалева // Вестник ЮУрГУ. Серия: Математическое моделирование и программирование. - 2015. - Т. 8, № 2. - С. 36-42.
13. Садин, Д.В. Применение схемы с настраиваемыми диссипативными свойствами к расчету течений газа с развитием неустойчивости на контактной границе / Д.В. Садин // Научно-технический вестник информационных технологий, механики и оптики. -2018. - Т. 18, № 1. - С. 153-157.
14. Gottlieb, S. Total Variation Diminishing Runge-Kutta Schemes / S. Gottlieb, C.W. Shu // Mathematics of Computation. - 1998. - V. 67, № 221. - P. 73-85.
15. Садин, Д.В. Сравнение разностной схемы с настраиваемыми диссипативными свойствами и схемы WENO на примере одномерных задач динамики газа и газовзвесей / Д.В. Садин, С.А. Одоев // Научно-технический вестник информационных технологий, механики и оптики. - 2017. - Т. 17, № 4. - С. 719-724.
16. Иванов, А.С. Нестационарное истечение двухфазной дисперсной среды из цилиндрического канала конечных размеров в атмосферу / А.С. Иванов, В.В. Козлов, Д.В Садин // Известия РАН. Механика жидкости и газа. - 1996. - № 3. - С. 60-66.
17. Садин, Д.В. Решение жестких задач течений двухфазных сред со сложной волновой структурой / Д.В. Садин // Физико-химическая кинетика в газовой динамике. - 2014.
- Т. 15, № 4. - С. 1-17.
18. Haas, J.-F. Interaction of Weak Shock Waves with Cylindrical and Spherical Gas Inhomogeneities / J.-F. Haas, B. Sturtevant // Journal of Fluid Mechanics. - 1987. - V. 181.
- P. 41-76.
19. Shyue, K.-M. An Eulerian Interface Sharpening Algorithm for Compressible Two-Phase Flow: the Algebraic THINC Approach / K.-M. Shyue, F. Xiao // Journal of Computational Physics.
- 2014. - V. 268. - P. 326-354.
20. Quirk, J.J. On the Dynamics of a Shock-Bubble Interaction / J.J. Quirk, S. Karni // Journal of Fluid Mechanics. - 1996. - V. 318. - P. 129-163.
Дмитрий Викторович Садин, доктор технических наук, профессор, Военно-космическая академия им. А.Ф. Можайского (г. Санкт-Петербург, Российская Федерация), sadin@yandex.ru.
Поступила в редакцию 27 сентября 2018 г.
MSC 76N15 DOI: 10.14529/mmp190209
A MODIFICATION OF THE LARGE-PARTICLE METHOD TO A SCHEME HAVING THE SECOND ORDER OF ACCURACY IN SPACE AND TIME FOR SHOCKWAVE FLOWS IN A GAS SUSPENSION
D. V. Sadin, Mozhaisky Military Space Academy, Saint-Petersburg, Russian Federation, sadin@yandex. ru
We develop the previously proposed approach of constructing difference schemes for solving stiff problems of shockwave flows of heterogeneous media using an implicit noniterative algorithm for calculating interactions between the phases. The large particle method is modified to a scheme having the second order of accuracy in time and space on smooth solutions. At the first stage, we use the central differences with artificial viscosity of TVD type. At the second stage, we implement TVD-reconstruction by weighted linear combination of upwind and central approximations with flow limiters. The scheme is supplemented by a two-step Runge-Kutta method in time. The scheme is K-stable, i.e. the time step does not depend on the intensity of interactions between the phases, but is determined by the Courant number for a homogeneous system of equations (without source terms). We use test problems to confirm the monotonicity, low dissipation, high stability of the scheme and convergence of numerical results to the exact self-similar equilibrium solutions in a gas suspension. Also, we show the scheme capability for numerical simulation of physical instability and turbulence. The method can be used for flows of gas suspensions having complex structure.
Keywords: gas-suspension flow; stiff problem; difference scheme; stability; accuracy.
References
1. Belotserkovskii O.M., Davydov Yu.M. A Non-Stationary "Coarse Particle" Method for Gas-Dynamical Computations. USSR Computational Mathematics and Mathematical Physics, 1971, vol. 11, no. 1, pp. 241-271. DOI: 10.1016/0041-5553(71)90112-1
2. 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. 16, no. 7, pp. 180-192. DOI: 10.1016/0041-5553(77)90183-5
3. Nigmatulin R.I. Dinamika mnogofaznyh sred [Dynamics of Multiphase Media]. Moscow, Nauka, 1987. (in Russian)
4. Sadin D.V., A Modified Large-Particle Method for Calculating Unsteady Gas Flows in a Porous Medium. Computational Mathematics and Mathematical Physics, 1996, vol. 36, no. 10, pp. 1453-1458.
5. Sadin D.V., A Method for Computing Heterogeneous Wave Flows with Intense Phase Interaction. Computational Mathematics and Mathematical Physics, 1998, vol. 38, no. 6, pp. 987-993.
6. Sadin D.V. 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, 2002, vol. 43, no. 2, pp. 286-290. DOI: 10.1023/A:1014714012032
7. Sadin D.V. TVD Scheme for Stiff Problems of Wave Dynamics of Heterogeneous Media of Nonhyperbolic Nonconservative Type. Computational Mathematics and Mathematical Physics, 2016, vol. 56, no. 12, pp. 2068-2078. DOI: 10.1134/S0965542516120137
8. Sadin D.V. Schemes with Customizable Dissipative Properties as Applied to Gas-Suspensions Flow Simulation]. Mathematical Models and Computer Simulations, 2017, vol. 29, no 12, pp. 89-104. (in Russian)
9. Saurel R., Abgrall R. A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows. Journal of Computational Physics, 1999, vol. 150, no. 2, pp. 425-467. DOI: 10.1006/jcph.1999.6187
10. Saurel R., Le Metayer O., Massoni J., Gavrilyuk S. Shock Jump Relations for Multiphase Mixtures with Stiff Mechanical Relaxation. Shock Waves, 2007, vol. 16, no. 3, pp. 209-232. DOI: 10.1007/ s00193-006-0065-7
11. Saurel R., Petitpas F., Berry R.A. Simple and Efficient Relaxation Methods for Interfaces Separating Compressible Fluids, Cavitating Flows and Shocks in Multiphase Mixtures. Journal of Computational Physics, 2009, vol. 228, no. 5, pp. 1678-1712. DOI: 10.1016/j.jcp.2008.11.002
12. Grishchenko D.S., Kovalev Yu.M., Kovaleva E.A. Modifcation 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)
13. Sadin D.V. 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, 2018, vol. 18, no. 1, pp. 153-157. DOI: 10.17586/2226-1494-2018-18-1-153-157 (in Russian)
14. Gottlieb S., Shu C.-W. Total Variation Diminishing Runge-Kutta Schemes. Mathematics of Computation, 1998, vol. 67, no. 221, pp. 73-85. DOI: 10.1090/S0025-5718-98-00913-2
15. Sadin D.V., Odoev S.A. Comparison of Difference Scheme with Customizable Dissipative Properties and WENO Scheme in the Case of One-Dimensional Gas and Gas-Particle Dynamics Problems. Scientific and Technical Journal of Information Technologies, Mechanics and Optics, 2017, vol. 17, no. 4, pp. 719-724. DOI: 10.17586/2226-1494-2017-17-4-719-724 (in Russian)
16. Ivanov A.S., Kozlov V.V., Sadin D.V. Unsteady Flow of a Two-Phase Disperse Medium from a Cylindrical Channel of Finite Dimensions Into the Atmosphere. Fluid Dynamics, 1996, vol. 31, no. 3, pp. 386-391. DOI: 10.1007/BF02030221
17. Sadin D.V. Stiff Problems of a Two-Phase Flow with a Complex Wave Structure. Physical-Chemical Kinetics in Gas Dynamics, 2014, vol. 15, no. 4, pp. 1-17. (in Russian)
18. Haas J.-F., Sturtevant B. Interaction of Weak Shock Waves with Cylindrical and Spherical Gas Inhomogeneities. Journal of Fluid Mechanics, 1987, vol. 181, pp. 41-76. DOI: 10.1017/S0022112087002003
19. Shyue K.-M., Xiao F. An Eulerian Interface Sharpening Algorithm for Compressible Two-Phase Flow: the Algebraic THINC Approach. Journal of Computational Physics, 2014, vol. 268, pp. 326-354. DOI: 10.1016/j.jcp.2014.03.010
20. Quirk J.J., Karni S. On the Dynamics of a Shock-Bubble Interaction. Journal of Fluid Mechanics, 1996, vol. 318, pp. 129-163. DOI: 10.1017/S0022112096007069
Received September 27, 2018
