Application of computational algorithms with higher order of accuracy to the modeling of two-dimensional problems on development of hydrodynamic instability

1.      Введение

Неустойчивость Рихтмайера-Мешкова (далее $–$ НРМ) возникает при прохождении ударной волны через разделительную границу между жидкостями разной плотности. Она инициализирована накоплением завихренности на границе раздела из-за сдвига градиентов давления и плотности на границе ударной волны и материала.

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

Исследования, посвященные численному моделированию развития неустойчивости Рихтмайера-Мешкова на основе уравнений Эйлера, не учитывали влияние газового взаимопроникновения (например, работы [1-2]). Кроме того, было отмечено, что замена ступенчатого профиля плотности на контактном разрыве непрерывным распределением в слое конечной ширины может снизить скорость роста возмущений на начальной стадии развития неустойчивости Рихтмайера-Мешкова. Этот факт подтвержден в исследованиях, проведенных в работах [3-4]. Также, в работе [5] было указано на необходимость использования моделей многокомпонентных смесей для описания разрушения контактной границы и образования области смеси.

Современные стандарты качества математического моделирования гидродинамических неустойчивостей требуют проведения расчетов на подробных сетках (до нескольких миллионов ячеек и более) с использованием высокоточных схем. Такие схемы, как правило, основанные на методах высокого порядка аппроксимации, активно разрабатываются и исследуются в настоящее время [6].

Актуальные и эффективные схемы для решения задач о развитии неустойчивости Рихтмайера-Мешкова включают в себя различные методы численного моделирования и анализа таких процессов. Некоторые из них включают в себя: метод конечных элементов [7], метод конечных разностей [8-9], метод сглаженных частиц [10] и т.д.

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

В данной работе исследовались две задачи о набегании ударной волны на область из более тяжелого газа с использованием схем второго и пятого порядка точности.

2.      Математическая модель

Будем рассматривать двумерную систему уравнений двухкомпонентной газовой динамики, записанную в консервативной форме:

                                                              $\frac{\partial U}{\partial t}+\frac{\partial F\mathrm{<mo\; stretchy="false">(</mo>}U\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x}+\frac{\partial G\mathrm{<mo\; stretchy="false">(</mo>}U\mathrm{<mo\; stretchy="false">)</mo>}}{\partial y}\mathrm{<mo>=</mo>0,}$ (2.1)

 где

                                    $U\mathrm{<mo>=</mo>}\left(\begin{array}{c}\rho \\ \rho u\\ \rho v\\ E\\ \rho Y\end{array}\right)\mathrm{,}F\mathrm{<mo\; stretchy="false">(</mo>}U\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\left(\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uv\\ \mathrm{<mo\; stretchy="false">(</mo>}E+p\mathrm{<mo\; stretchy="false">)</mo>}u\\ \rho uY\end{array}\right)\mathrm{,}G\mathrm{<mo\; stretchy="false">(</mo>}U\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\left(\begin{array}{c}\rho v\\ \rho uv\\ \rho v^2+p\\ \mathrm{<mo\; stretchy="false">(</mo>}E+p\mathrm{<mo\; stretchy="false">)</mo>}v\\ \rho vY\end{array}\right)\mathrm{.}$

 Здесь $\rho$ $–$ плотность жидкости, $v\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}u\mathrm{,}v\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ вектор скорости, $p$ $–$ давление и $E\mathrm{<mi>\; </mi>}\text{}\mathrm{<mo>=</mo><mi>\; </mi>}\text{}\rho \mathrm{<mo\; stretchy="false">(</mo>}\epsilon \mathrm{<mi>\; </mi>}+\mathrm{<mi>\; </mi>}\frac{u^2+v^2}{2}\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ полная энергия, $\epsilon$ $–$ удельная внутренняя энергия идеального газа.

Система (2.1) замыкается уравнением состояния $p\mathrm{<mo>/</mo>}\rho \mathrm{<mo>=</mo>}RT\mathrm{<mo>/</mo>}M\mathrm{,}R\mathrm{<mo>/</mo>}M\mathrm{<mo>=</mo>}{C}_p-C_v$, где $R$ $–$ универсальная газовая постоянная, $C_p\mathrm{,}{C}_v$ $–$ удельные теплоемкости смеси при постоянном давлении и постоянном объеме соответственно, $M$ $–$ молекулярная масса смеси. $C_p$ и $M$ вычисляются следующим образом:

                                                                  $C_p\mathrm{<mo>=</mo>}Y{C}_{p1}+\mathrm{<mo\; stretchy="false">(</mo>1}-Y\mathrm{<mo\; stretchy="false">)</mo>}{C}_{p2}\mathrm{,}$ (2.2)

                                                                         $\frac{1}{M}\mathrm{<mo>=</mo>}\frac{Y}{M_1}+\frac{1-Y}{M_2}\mathrm{,}$ (2.3)

где $Y$ $–$ концентрация, $M_i\mathrm{,}{C}_{pi}$ $–$ молекулярная масса и теплоемкость $i$ -й компоненты смеси соответственно ( $i\mathrm{<mo>=</mo>1,2}$ ).

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

3.      Вычислительный алгоритм

Построим дискретную модель для расчетной области прямоугольной формы $\Omega \mathrm{<mi>\; </mi>}\text{}\text{}\mathrm{<mo>=</mo><mi>\; </mi>}\text{}\text{}\mathrm{<mo\; stretchy="false">[</mo>0,}{L}_x\mathrm{<mo\; stretchy="false">]</mo><mi>\; </mi>}\text{}\text{}\text{}\times \mathrm{<mi>\; </mi>}\text{}\text{}\text{}\mathrm{<mo\; stretchy="false">[</mo>0,}{L}_y\mathrm{<mo\; stretchy="false">]</mo>}$. Для этого область $\Omega$ заменим ортогональной сеткой, равномерной по каждому направлению:

                                                                          ${\omega }_{\Delta }\mathrm{<mo>=</mo>}{\omega }_{{\Delta }_x}\times {\omega }_{{\Delta }_y}\mathrm{,}$

 где

                      $\begin{array}{l}{\omega }_{{\Delta }_x}\mathrm{<mo>=</mo><mo>\{</mo>}{\Delta }_i\mathrm{,}i\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{N}_x\mathrm{,}{\Delta }_i\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>}{x}_{i-\mathrm{1<mo>/</mo>2}}\mathrm{,}{x}_{i+\mathrm{1<mo>/</mo>2}}\mathrm{<mo\; stretchy="false">]</mo>,<mo>|</mo>}{\Delta }_i\mathrm{<mo>|</mo><mo>=</mo>}{h}_x\mathrm{,}{h}_x{N}_x\mathrm{<mo>=</mo>}{L}_x\mathrm{<mo>\}</mo>,}\hfill \\ {\omega }_{{\Delta }_y}\mathrm{<mo>=</mo><mo>\{</mo>}{\Delta }_j\mathrm{,}j\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{N}_y\mathrm{,}{\Delta }_j\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>}{y}_{j-\mathrm{1<mo>/</mo>2}}\mathrm{,}{y}_{j+\mathrm{1<mo>/</mo>2}}\mathrm{<mo\; stretchy="false">]</mo>,<mo>|</mo>}{\Delta }_j\mathrm{<mo>|</mo><mo>=</mo>}{h}_y\mathrm{,}{h}_y{N}_y\mathrm{<mo>=</mo>}{L}_y\mathrm{<mo>\}</mo>.}\hfill \end{array}$ (3.1)

Для аппроксимации системы (2.1) используем нелинейную консервативную дифференциально-разностную схему повышенного порядка точности:

                                               $\frac{d{U}_{ij}}{dt}+\frac{F_{i+\mathrm{1<mo>/</mo>2}j}-F_{i-\mathrm{1<mo>/</mo>2}j}}{h_x}+\frac{G_{ij+\mathrm{1<mo>/</mo>2}}-G_{ij-\mathrm{1<mo>/</mo>2}}}{h_y}\mathrm{<mo>=</mo>0,}$ (3.2)

 где $U_{ij}\mathrm{<mo>=</mo>}\frac{1}{h_x{h}_y}{\displaystyle \underset{x_{i-\mathrm{1<mo>/</mo>2}}}{\iint }}{\displaystyle \underset{y_{j-\mathrm{1<mo>/</mo>2}}}{\overset{y_{j+\mathrm{1<mo>/</mo>2}}}{\int }}}Udxdy$ $–$ соотнесенные к центру ячейки усредненные значения консервативных газодинамических переменных, $F_{i+\mathrm{1<mo>/</mo>2}j}\mathrm{,}{G}_{ij+\mathrm{1<mo>/</mo>2}}$ $–$ дискретные потоки на соответствующих границах между ячейками, являющиеся функциями двух переменных

                                                                $F_{i+\mathrm{1<mo>/</mo>2}j}\mathrm{<mo>=</mo>}\widehat{F}\mathrm{<mo\; stretchy="false">(</mo>}{U}_{i+\mathrm{1<mo>/</mo>2}j}^L\mathrm{,}{U}_{i+\mathrm{1<mo>/</mo>2}j}^R\mathrm{<mo\; stretchy="false">)</mo>,}$ (3.3)

                                                                $G_{i+\mathrm{1<mo>/</mo>2}j}\mathrm{<mo>=</mo>}\widehat{G}\mathrm{<mo\; stretchy="false">(</mo>}{U}_{ij+\mathrm{1<mo>/</mo>2}}^L\mathrm{,}{U}_{ij+\mathrm{1<mo>/</mo>2}}^R\mathrm{<mo\; stretchy="false">)</mo>,}$ (3.4)

 для которых выполнено условие согласования:

                                                                      $\widehat{F}\mathrm{<mo\; stretchy="false">(</mo>}{U}_{ij}\mathrm{,}{U}_{ij}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}F\mathrm{<mo\; stretchy="false">(</mo>}{U}_{ij}\mathrm{<mo\; stretchy="false">)</mo>,}$ (3.5)

                                                                      $\widehat{G}\mathrm{<mo\; stretchy="false">(</mo>}{U}_{ij}\mathrm{,}{U}_{ij}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}G\mathrm{<mo\; stretchy="false">(</mo>}{U}_{ij}\mathrm{<mo\; stretchy="false">)</mo>.}$ (3.6)

 Здесь $U_{i+\mathrm{1<mo>/</mo>2}j}^L\mathrm{,}{U}_{i+\mathrm{1<mo>/</mo>2}j}^R$ $–$ «левые» и «правые» значения вектора $U$ на грани между $i$ и $i+1$ ячейками, для которой вычисляется поток $F_{i+\mathrm{1<mo>/</mo>2}j}$, $m\mathrm{<mo>=</mo>1,2,3}$; $U_{ij+\mathrm{1<mo>/</mo>2}}^L\mathrm{,}{U}_{ij+\mathrm{1<mo>/</mo>2}}^R$ $–$ «левые» и «правые» значения вектора $U$ на грани между $j$ и $j+1$ ячейками. Для того, чтобы вычислить значения вектора $U$ на указанных гранях между ячейками введем новый вектор переменных $Q\mathrm{<mo>=</mo>}Q\mathrm{<mo\; stretchy="false">(</mo>}U\mathrm{<mo\; stretchy="false">)</mo>}$. Проведем его интерполяцию на грань между ячейками и затем пересчитаем искомое значение вектора $U\mathrm{<mo>=</mo>}U\mathrm{<mo\; stretchy="false">(</mo>}Q\mathrm{<mo\; stretchy="false">)</mo>}$ на данной грани. В расчетах полагалось $Q\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}\rho \mathrm{,}p\mathrm{,}u\mathrm{,}v\mathrm{,}Y\mathrm{<mo\; stretchy="false">)</mo>}$, дискретные потоки вычислялись по схемам Лакса-Фридрихса (LF) [11] и Хартена-Лакса-ван Лира с учетом контактного разрыва (HLLC) [12].

Для интерполяции значений $Q$ на грани между ячейками использовались схема TVD с лимитером $minmod$ [9] (далее $–$ TVD2) и схема WENO5 [13].

Дискретизация по времени проводилась с использованием TVD-схемы Рунге-Кутта 3-го порядка [13]. А именно, для уравнения вида

                                                                             $\frac{\partial U}{\partial t}\mathrm{<mo>=</mo>}{L}_h\mathrm{<mo\; stretchy="false">(</mo>}U\mathrm{<mo\; stretchy="false">)</mo>,}$ (3.7)

 где $L_h\mathrm{<mo\; stretchy="false">(</mo>}U\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ пространственный разностный оператор из (3.2), используем:

                 $U^{\mathrm{<mo>*</mo>}}\mathrm{<mo>=</mo>}{U}^n+\Delta t\cdot L_h\mathrm{<mo\; stretchy="false">(</mo>}{U}^n\mathrm{<mo\; stretchy="false">)</mo>,}$ (3.8)

                 $U^{\mathrm{<mo>*</mo><mo>*</mo>}}\mathrm{<mo>=</mo>}\frac{3}{4}{U}^n+\frac{1}{4}{U}^{\mathrm{<mo>*</mo>}}+\frac{1}{4}\Delta t\cdot L_h\mathrm{<mo\; stretchy="false">(</mo>}{U}^{\mathrm{<mo>*</mo>}}\mathrm{<mo\; stretchy="false">)</mo>,}$ (3.9)

                 $U^{n+1}\mathrm{<mo>=</mo>}\frac{1}{3}{U}^n+\frac{2}{3}{U}^{\mathrm{<mo>*</mo><mo>*</mo>}}+\frac{2}{3}\Delta t\cdot L_h\mathrm{<mo\; stretchy="false">(</mo>}{U}^{\mathrm{<mo>*</mo><mo>*</mo>}}\mathrm{<mo\; stretchy="false">)</mo>.}$ (3.10)

4.      Постановка задач

Задача 1.

 Рассматривается задача о развитии неустойчивости Рихтмайера-Мешкова при прохождении ударной волны через прямоугольную неоднородность из тяжелого газа. Начально-краевые условия взяты согласно эксперименту [14] и расчетам, описанным в работе [15], в которой используется схема КАБАРЕ.

Рис. 4.1. Схема расчетной области для задачи 1

Fig 4.1. Scheme of the computational domain for problem 1 

Расчетная область $–$ прямоугольная $\mathrm{<mo\; stretchy="false">[</mo>0,}{L}_x\mathrm{<mo\; stretchy="false">]</mo>}\times \mathrm{<mo\; stretchy="false">[</mo>0,}{L}_y\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo><mo\; stretchy="false">[</mo>0,0.45<mo\; stretchy="false">]</mo>}\times \mathrm{<mo\; stretchy="false">[</mo>0,0.2<mo\; stretchy="false">]</mo>}$, изображена на Рис. 4.1. В начальный момент времени в примыкающей к стенке прямоугольной подобласти $\mathrm{<mo\; stretchy="false">[</mo>0.1,0.25<mo\; stretchy="false">]</mo>}\times \mathrm{<mo\; stretchy="false">[</mo>0,0.1<mo\; stretchy="false">]</mo>}$ (подобласть II на Рис. 4.1) находится покоящийся тяжелый газ $–$ фторид серы VI ( $S{F}_6$, элегаз), в остальной части находится покоящийся воздух (подобласть I на Рис. 4.1), оба газа находятся в статическом равновесии. На левой грани задается условие входа ударной волны, на остальных гранях $–$ адиабатические стенки с проскальзыванием. Параметры сред, принятые для вычислительного эксперимента, указаны в таблице 4.1.

Таблица 4.1. Данные для численного эксперимента для задачи 1

Table 4.1. Data for numerical experiment for problem 1 

Задача 2. Была принята следующая постановка задачи в соответствии с экспериментом, описанным в работе [16]. Расчетная область $–$ прямоугольная      $\mathrm{<mo\; stretchy="false">[</mo>0,}{L}_x\mathrm{<mo\; stretchy="false">]</mo>}\times \mathrm{<mo\; stretchy="false">[</mo>0,}{L}_y\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo><mo>=</mo><mo\; stretchy="false">[</mo>0,0.2<mo\; stretchy="false">]</mo>}\times \mathrm{<mo\; stretchy="false">[</mo>0,0.14<mo\; stretchy="false">]</mo>}$ (Рис. 4.2).

Рис. 4.2. Схема расчетной области для задачи 2

Fig 4.2. Scheme of the computational domain for problem 2

В начальный момент времени в квадратной подобласти II со стороной $a\mathrm{<mo>=</mo>0.0566}$ находится покоящийся тяжелый газ $–$ фторид серы VI ( $S{F}_6$, элегаз) в остальной части находится азот (подобласти I и III на Рис. 4.2). Фронт ударной волны, движущейся вправо с числом Маха 1.17, находится на расстоянии $L_{sw}\mathrm{<mo>=</mo>0.02}$ от левой границы. Область с тяжелым газом располагается перед фронтом ударной волны на расстоянии, соответствующем достижению ударной волной контактной границы за 4 мкс и по центру относительно оси $Oy$. Давление полагалось равным 101 325 Па, а температура $–$ равной 298 К. На левой грани задается условие втекания газа с параметрами за фронтом ударной волны, на остальных гранях $–$ условия свободного вытекания. Параметры газов приведены в таблице 4.2.

Таблица 4.2. Данные для численного эксперимента для задачи 2

Table 4.2. Data for numerical experiment for problem 2 

5.      Результаты расчетов

5.1. Задача 1 

На Рис. 5.1 показаны результаты расчетов различными методами на сетках размерности $360\times 160$, $720\times 320$ и $1080\times 480$. Видно, что использование солвера HLLC для вычисления дискретных потоков позволяет воспроизвести корректную картину течения даже на грубой сетке, а использование реконструкции WENO5 позволяет смоделировать достаточно подробную структуру течения.

В работе [14] приводятся результаты измерения положения границ области с элегазом в различные моменты времени. На основе этих данных было вычислено значение ширины зоны с элегазом в эксперименте и выполнено сравнение с результатами расчетов в данной работе. На Рис. 5.2 приводится сравнение изменения ширины зоны с элегазом для различных методов на последовательности измельчающихся сеток. Видно, что до момента ( $\approx 1.5$ мс) прохождения вторичной ударной волны, отраженной от правой границы, все методы демонстрируют приемлемое совпадение изменения ширины наблюдаемой зоны с элегазом. Далее во всех расчетах значение ширины занижено, за исключением результатов, полученных с использованием реконструкции WENO5, которые попадают в доверительный интервал экспериментальных данных.

Рис. 5.1. Задача 1 – картина течения на момент времени 4 046 мкс: a) эксперимент [14]; b) результаты расчетов (концентрация)

Fig 5.1. Problem 1  flow pattern at time 4 046  s: a) experiment [14]; b) calculation results (concentration)

Рис. 5.2. Задача 1 – изменение ширины области с элегазом: a) схема первого порядка точности, поток LF; b) схема первого порядка точности, поток HLLC; c) реконструкция TVD2, поток LF; d) реконструкция TVD2, поток HLLC; e) реконструкция WENO5, поток LF; f) реконструкция WENO5, поток HLLC

Fig 5.2. Task 1  changing the width of the area with SF6 gas: a) first-order accuracy circuit, LF flux; b) first order accuracy circuit, HLLC flux; c) reconstruction of TVD2, LF flux; d) TVD2 reconstruction, HLLC flux; e) WENO5 reconstruction, LF flux; f) WENO5 reconstruction, HLLC flux

 5.2 Задача 2

 На Рис. 5.3 показаны численные шлирен-картины по истечении 887 мкс с момента достижения ударной волной левой границы области II.

Рис. 5.3. Задача 2 – численные шлирен-картины на момент времени 877 мкс: a) эксперимент [16]; b) результаты расчетов

Fig 5.3. Problem 2 $–$ numerical schlieren pictures at time 877 $\mu$ s: $a$) experiment [16]; $b$) calculation results

Легко заметить, что с повышением порядка точности схемы увеличивается детализация картины течения. Схема на основе WENO позволяет разрешить более «тонкие» детали развития вихревых структур на границе раздела двух газов. Картины течения, полученные с помощью WENO-схем, имеют более закрученные вихревые структуры, чем в эксперименте [16] (Рис. 5.3 $a$ ). Это можно объяснить тем, что в модели не учитывается вязкость, которая увеличивала бы диссипацию вихревых структур.

На Рис. 5.4 представлена динамика изменения ширины области с тяжелым газом в сравнении с экспериментом.

Рис. 5.4. Изменение ширины области с тяжелым газом: a) сетка ; b) сетка ; c) сетка ; d) схема WENO5 с потоками LF и HLLC на сгущающихся сетках

Fig 5.4. Dynamics of the width of the region with heavy gas: $a$ ) grid $250\times 90$; $b$ ) grid $500\times 180$; $c$ ) grid $1000\times 360$; $d$ ) WENO5 scheme with LF and HLLC fluxes on condensed grids

Видно, что наиболее близко к эксперименту изменение ширины области с элегазом воспроизводит схема с дискретным потоком HLLC и реконструкцией WENO 5-го порядка точности.

6.      Заключение

В работе решены две задачи о развитии гидродинамической неустойчивости при набегании ударной волны на область с более плотным газом. Получены картины течения и проанализирована динамика изменения области с плотным газом. Расчеты проведены с использованием вычислительных алгоритмов на основе интегро-интерполяционного метода первого порядка точности и с реконструкцией решения на границах ячеек по схеме TVD второго порядка точности и по схеме WENO птого порядка точности. Для вычисления дискретных потоков на границах ячеек использовались потоки Лакса-Фридрихса и HLLC. Что наиболее близкие к эксперименту результаты демонстрирует схема с реконструкцией WENO и дискретными потоками HLLC.

Благодарности. Постановка расчетных задач, обработка и интерпретация результатов выполнена Жалниным Р. В. за счет средств РНФ (проект $№$ 23-11-00142).