Distinctive features of calculating interaction between shock wave and gas bubble in a finely dispersed gas suspension.pdf Потоки газовзвесей находят широкое применение в практических приложениях [1-3]. Для таких явлений характерно наличие существенного диапазона изменения газодинамических параметров течений, границы раздела сред, неоднородностей по плотности и крупных включений. Несомненный интерес представляют структурно-сложные задачи, которые отличают рефракция ударных волн на границе раздела сред и одновременно развитие неустойчивости и турбулентности. При численном моделировании ударно-волновых процессов в газовзвесях с мелкодисперсной фракцией возникает ряд проблем. Для таких течений характерно явление жесткости, связанное с существенным различием времен релаксации фаз. Применение традиционных явных алгоритмов расчета становится практически не пригодным ввиду неоправданно малого шага по времени. Для преодоления этой трудности более 20 лет назад предложены неявные и полунеявные методы расчета источниковых членов (межфазных взаимодействий) [4-6]. На модельных односкоростных и двухскоростных уравнениях в частных производных показано, что такой подход справедлив для широкого класса разностных схем [4, 5]. Продемонстрирована высокая устойчивость схем даже для некорректно поставленных разностных задач - возможность устойчивого проведения расчетов с числом Куранта больше единицы CFL>1. В частности, для линейных задач построены аппроксимирующие схемы, доказана их устойчивость при CFL > 1, но из этих двух свойств не следует сходимость [6]. Проведенный анализ схем для двухскоростной динамики гетерогенных сред показал, что для практически важных случаев достаточно неявно рассчитывать только быстрые компоненты решения (газовую фазу). Изучены влияние на допустимый шаг по времени (число Куранта) шага расчетной сетки, размеров частиц газовзвеси, концентрации дисперсной фазы и других величин [5]. Введено понятие K-устойчивости схемы - условия устойчивости не зависят от интенсивности межфазных взаимодействий и размеров сетки [7, 8]. Другой подход основан на расщеплении законов сохранения газовзвесей на однородную систему уравнений (без источников) и введение корректирующего шага для релаксационных процессов в виде обыкновенных дифференциальных уравнений [9]. Однако в работе [10] показано, что указанная стратегия построения разностной схемы может приводить к существенному дисбалансу градиентно-потоковых и источниковых членов уравнений, например, при получении численного решения методом установления. Интерес в исследованиях при численном моделировании ударно-волновых течений газовзвесей вызывает явление жесткости [11]. Весьма важным требованием, предъявляемым к разностным схемам, является возможность разрешения деталей ударной волновой картины и часто сопровождаемых явлений развития неустойчивости на контактных поверхностях и турбулентности. Схемы первого порядка точности являются избыточно диссипативными и их применение требует значительных временных затрат для получения требуемого качества решения для структурно-сложных течений газа и газовзвесей. Для сохранения монотонности при повышении порядка аппроксимации схемы используют схемы высокой разрешающей способности типа Годунова [9, 12], TVD [10], WENO-типа [13], с адаптивной искусственной вязкостью [14] и др. Еще ряд проблем связан с численным моделированием взаимодействия волн с границами раздела сред с различными уравнениями состояния и сходимостью численных решений к равновесным течениям при уменьшении размеров частиц газовзвеси. Газодинамические течения в рамках сплошной среды в общем случае описываются уравнениями Навье-Стокса. Численное моделирование таких течений является корректным при условии - численная диффузия существенно меньше физической. В настоящей работе используется подход [13, 15] в рамках уравнения Эйлера, например для оценки разрешения сеток, обеспечивающих в последующих исследованиях корректное моделирование на основе уравнений Навье-Стокса при заданном числе Рейнольдса [13]. Постановка задачи Для описания ударно-волновых процессов в газовзвесях запишем уравнения двухскоростной двухтемпературной динамики двухфазной среды [16] при известных допущениях в виде [14] ∂q + V d G + B (V d F ) = H (q ), q=[ρ1,ρ2,ρ1v1,ρ2v2,ρ2e2,ρ1E1+ρ2K2]T; (1) G=[ρ1v1,ρ2v2,ρ1v1v1,ρ2v2v2,ρ2e2v2,ρ1E1v1+ρ2K2v2]T; F = [0,0, р, р, 0, p (α1Γ1 + α2Г2 )]t , H = [0,0, -Fμ, Fμ, Qt , -Qt ]Т ; Vd = diag (V∙, V∙,V,V,V∙,V∙) , B = diag [1,1, α1, а 2,1,1]; Pi = Pi αi (i = 1,2), E1 = e1 + v12 /2, K2 = v2/2 . Здесь и далее индексы 1 и 2 внизу относятся соответственно к параметрам несущей и дисперсной фаз, индекс ° сверху относится к истинным значениям плотности; V - оператор Гамильтона. Через αi, ρi, v1, E1, e1, р обозначены объемная доля, приведенная плотность, вектор скорости, полная и внутренняя энергии единицы массы i-й фазы, давление газа; Fμ, Qt - соответственно вязкая составляющая силы межфазного взаимодействия, мощность теплообмена между газом и частицами в единице объема; t - время. Особенности расчета взаимодействия ударной волны с газовым пузырем 101 Для замыкания системы (1) используем уравнения состояния идеального кало-рически совершенного газа и несжимаемых твердых частиц: p = (γ1 - 1)ρ1°e1, e1 = cvT1, e^ = c2T2, {γ1,cv,c2^,ρ2! ≡ const, где T1, T2 - температура несущей фазы и частиц; γ1,cv - показатель адиабаты и удельная теплоемкость газа при постоянном объеме; c2 - удельная теплоемкость частиц. Интенсивности межфазного трения и теплообмена Fμ,QT задаются на основе экспериментальных соотношений [17-19]. Начальные условия задаются применительно к рассматриваемым ниже задачам. Граничные условия: на стенках - отражения, на внешних границах (во всех рассматриваемых ниже задачах слева и справа расчетной области) - экстраполяция. Метод решения Для численного решения используем схему с настраиваемыми диссипативными свойствами [14, 20] в конечно-объемной реализации [21] с аппроксимацией второго порядка по времени и пространству на гладких решениях - CDP2 (customizable dissipative properties). Переход на следующий временной слой выполняется двухшаговым TVD-методом Рунге-Кутты [22]: q(1) = qk +τL(qk) , qk+1 = 0.5(qk + q‘-k)) + 0.5τL(q(1)) . При этом на каждом шаге производится расщепление пространственного разностного оператора L (q) по физическим процессам с использованием искусственной вязкости Qnk±1/2 типа Христенсена [23] с ограничителями вязкости ψv [20, 21]: - δH(qn ■ )т = qn +(1 -δ)H(q^ )т-Вк„ (FF^h,, - irk_12 )т/h , qk+1/2 qn 1) где F'n±1/2 =[°,°, pn±1∕2, pn±1∕2,0, pn±1∕2 (а1^1 +а2V2 )n±1/2 ] ∙> pn±1∕2 = pn±1∕2 + 6n±1∕2 ; = qk+1/2-τ[Vd G ]T+;^. qn(1) 2) Здесь сится полуцелые n±1∕2 к ее граням; [VdG]k+1/2 n т, h - временные и пространственные шаги сетки; верхний индекс k отно-к временному слою; нижний целый индекс n относится к центру ячейки, а - TVD-реконструкция потоков путем взвешенной линейной комбинации противопоточной и центральной аппроксимаций конвективных членов с ограничителями потоков ψ f [21]; δ- весовой множитель источниковых членов (δ=1 неявный учет межфазных взаимодействий, δ=0 - явный). Высокая устойчивость обеспечивается учетом межфазных взаимодействий на верхнем временном слое (δ=1), что ранее было предложено и обосновано в работах [4, 5, 14 и др.]. Искомые функции рассчитываются явно (без итераций) при линейной зависимости источниковых членов от решения. Однако при использо- 102 Д.В. Садин, В.А. Давидчук вании эмпирических критериальных соотношений такие зависимости - нелинейные. Например, сила межфазного трения определяется суммой слагаемых, зависящих от разности скоростей фаз в различных степенях. Для исключения итерационных процедур используем линеаризацию источников, учитывая неявно линейную часть [5]. Как показал опыт расчетов течений газовзвесей в широком диапазоне параметров потоков и размеров частиц (от "00 нм и выше), запас устойчивости алгоритма при таком подходе сохраняется. Число Куранта для всех задач в настоящей работе принято CFL = 0.4, а ограничители искусственной вязкости и потоков - Van Leer [2"]. Верификация метода Особенностью рассматриваемой задачи взаимодействия ударной волны с газовым пузырем в мелкодисперсной газовзвеси является сложная ударно-волновая структура и развивающаяся неустойчивость (турбулентность) на контактном разрыве. С целью проверки возможности метода CDP2 для разрешения тонких деталей течения, диссипативных свойств и монотонности схемы рассмотрим тестовые задачи. Первая из них - задача распада разрыва в газовзвеси, имеющая точное автомодельное решение в равновесном приближении [24, "4]. Тестовая задача решается на интервале х / L ∈ (0,") с начальным разрывом в точке х0 / L = 0.48 и на отрезке времени t∈ [0, t" ] , где L" = " м - размер расчетной сетки с разрешением h/L" ="/500; t" =".2 мс - время окончания счета. Начальные условия задачи в системе СИ представлены в табл. ". Газодинамические величины помечены индексом L - слева от разрыва, а индексом R - справа. Двухфазная смесь в начальный момент находится в состоянии покоя. Частицы - монодисперсные, диаметром d = 0." мкм , поровый газ - воздух с показателем адиабаты γ = 7/5. После распада начального разрыва влево по газовзвеси распространяется волна разрежения r (рис. "), а вправо ударная волна s, разделенные комбинированным разрывом (скачком пористости) c и зоной постоянного течения между ними. Для заданных мелких частиц времена скоростной и тепловой релаксаций малы, много меньше шага по времени, поэтому течение практически равновесное дится к точному решению (рис. ", сплошные кривые). и схо- Таблица " Начальные условия тестовой задачи распада разрыва в газовзвеси α2L РЇі pl α2R ρ"R pr 0.00" "".89"9 "06 0.00" "."89"9 "05 Второй одномерный тест представляет собой модификацию известной задачи Карни и Кёрка [25]. В нашем случае ударная волна распространяется по газовзве-си и взаимодействует с ограниченной областью «чистого» газа (воздуха). Особенностью данной задачи является прохождение ударной волны и ее взаимодействие со средами, описываемыми различными уравнениями состояний. Кроме того, вычислительные трудности при решении этой задачи связаны с большим перепадом плотности. Особенности расчета взаимодействия ударной волны с газовым пузырем 103 0.75 0.50 p / Pl , ρ / Pl a Рис. 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 c 2 1 0.25 r s В начальный момент времени ударная волна, удовлетворяющая соотношениям Рэнкина - Гюгонио в равновесной газовзвеси [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 Начальные условия модифицированной задачи Карни и Кёрка Характерные области α2 РЇ 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-i0 1.18919 0 100000 104 Д.В. Садин, В.А. Давидчук Рис. 2. Сравнение численных решений (пунктирные кривые) с эталонными (сплошные линии) модифицированной задачи Карни и Кёрка в газовзвеси в момент времени t2 = 2 мс: a -плотность смеси, отнесенная к этой величине за ударной волной; b - скорость смеси, отнесенная к скорости звука в газе перед ударной волной Fig. 2. Comparison of numerical (dashed lines) and benchmark solutions (solid lines) to a modified Karni 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 Еще одна важная особенность связана с пространственными эффектами интерференции волн разрежения, ударных волн и контактных разрывов, разделяющих два газа с различными показателями адиабаты. Рассмотрим показательную с этой точки зрения задачу о взаимодействии ударной волны с цилиндрическим пузырем гелия, которая нашла широкое применение при тестировании разностных схем. Ударная волна с числом Маха 1.22 движется влево по каналу с поперечным размером 8.9 см и встречает пузырь диаметром 5 см смеси гелия и воздуха с массовой долей 28 %. Перед ударной волной воздух и воздушно-гелиевая смесь находятся в состоянии покоя при температуре 293.23 К и давлении 101325 Па. Результаты расчета по схеме CDP2 в виде численных шлирен-изображений функции градиента плотности газа показаны на рис. 3, a, b, c. Время отсчитывается от момента прихода ударной волны на границу пузыря. Для сравнения приведем численные данные, полученные по схеме WENO-типа четвертого порядка точности по пространству и третьего порядка по времени [15] (рис. 3, d, e, f). Расчеты выполнены при одинаковых исходных данных, числе Куранта и дискретизации сеток с шагом h=50 мкм. Численные решения хорошо согласуются между собой. Схема CDP2 демонстрирует возможности по разрешению тонких деталей ударно-волновой структуры, обладает малой диссипацией при сохранении монотонности решения. При этом вычислительные затраты существенно меньше, чем при использовании цитируемой WENO-схемы. Результаты расчета Продолжим рассмотрение задачи распада разрыва в газовзвеси при изменении размеров частиц, когда эффекты неравновесности становятся существенными. Задача решается при исходных данных (см. табл. 1), а диаметры частиц заданы 10 и 100 мкм. После распада разрыва возникает релаксирующее течение газовзвеси. Особенности расчета взаимодействия ударной волны с газовым пузырем 105 Рис. 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 2a(,2d( 62; 2^(,2e( 80; and 2с'(,2f( 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, штриховые кривые). Комбинированный разрыв c, который образовывался в равновесной газовзвеси (рис. 1, a), распадается на контактный разрыв в поровом газе i и скачок пористости c (рис. 4, a). В смеси с частицами d=10 мкм ударная волна вырождается в звуковую волну (рис. 4, штрихпунктирные кривые). При уменьшении размеров частиц зоны релаксации монотонно уменьшаются, а численное решение сходится к точному (рис. 4, сплошные линии). В практически важных случаях ударно-волновые возмущения распространяются в газовзвесях, неоднородных по плотности. Поэтому задача взаимодействия ударной волны с газовым пузырем представляет как теоретический интерес, так и прикладное значение. Геометрические параметры соответствуют задаче взаимодействия ударной волны с гелиевым пузырем. А начальные условия - см. табл. 2 (за исключением скорости, которая имеет обратный знак). Задача решалась на равномерной сетке 4500×534 ячеек с одинаковыми шагами по пространственным координатам. Плоская ударная волна распространяется в канале по мелкодисперсной (d = 0.1 мкм) газовзвеси и в некоторый момент времени встречает на своем пути неоднородность в виде воздушного пузыря. При этом возникает распад разрыва: с ударной волной внутри пузыря и ее рефракцией во внешнюю область, а также отраженной волной разрежения в противоположном направлении (см. рис. 5, а). Рис. 4. Сравнение численных решений (штрихпунктирные - d=10 мкм, пунктирные кривые - d = 100 мкм) с точными равновесными решениями (сплошные линии) задачи о распаде разрыва в газовзвеси в момент времени tf =2 мс: a - относительные давления 1 и плотности смеси 2; b - относительные скорости частиц 3, 5 (d=10 мкм) и газа 4, 6 (d = 100 мкм) Fig. 4. Comparison of numerical solutions (dotted-and-dashed lines at d=10 μm; dashed lines at d = 100 μm) and exact equilibrium solutions (solid lines) to the problem of decay of discontinuity in a gas suspension at the time instant of tf =2 ms: (a) the relative values of 1, pressure and 2, density of mixture; (b) the relative velocity of 3, 5, particles (d=10 μm) and 4, 6, gas (d = 100 μm) Рис. 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 μs Особенности расчета взаимодействия ударной волны с газовым пузырем 107 Скорость звука в газе выше, чем в газовзвеси. Поэтому ударная волна внутри пузыря движется быстрее, чем в окружающей газодисперсной среде. В момент времени от начала расчета t=250 мкс (рис. 5, b) взаимодействующая ударная волна прошла цилиндрическую неоднородность с образованием тройной точки. На рис. 5, c отмечается начало зарождающей неустойчивости Рихтмайера - Мешкова на поверхности комбинированного разрыва. Последующие моменты времени характеризуются более сложной ударно-волновой картиной и развитием крупномасштабной турбулентности, что хорошо видно на приведенных шлирен-изобра-жениях (рис. 5, d-f). Выводы Исследование проведено в рамках модели двухскоростной двухтемпературной динамики газовзвеси. Для расчета использована разностная схема со вторым порядком точности по пространству и времени CDP2, пригодная для решения жестких задач. Высокая устойчивость схемы обеспечена неявным учетом межфазных взаимодействий. Схема верифицирована на тестовой задаче распада разрыва в га-зовзвеси в сравнении с точным автомодельным решением в равновесной газовзве-си, модифицированной задаче Карни и Кёрка, а также сопоставлены с данными других авторов по взаимодействию ударной волны с цилиндрической неоднородностью в виде гелиевого пузыря. Тесты подтвердили малую диссипацию и монотонность численных решений, полученных по схеме CDP2. Рассмотрены релаксирующие течения газовзвесей при изменении размеров частиц и их сходимость к равновесным решениям. Численно решена задача о распространении ударной волны в мелкодисперсной газовзвеси и ее взаимодействии с пузырем воздуха. Рассмотрены ударно-волновая картина, развитие неустойчивости Рихтмайера -Мешкова на поверхности раздела сред и образование крупномасштабной турбулентности. Проведенные исследования позволяют оценить численную диссипацию на заданных разрешениях сеток для корректного моделирования течений га-зовзвесей в рамках Навье - Стокса при заданном числе Рейнольдса, что является направлением дальнейших исследований.

Скачать электронную версию публикации