Mathematical modeling of turbulent mixing in gas systems with a chevron contact boundary using NUT3D, BIC3D, EGAK, and MIMOSA numerical codes

全文:

1. ВВЕДЕНИЕ

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

На начальном этапе имеющиеся на контактной границе веществ малые возмущения под действием неустойчивостей начинают расти и довольно быстро выходят на насыщение и нелинейную стадию, на которой происходит взаимодействие и разрушение различных мод и формирование зоны турбулентного перемешивания. Законы развития этой зоны перемешивания, а также варианты ее взаимодействия с классическими сильными и слабыми разрывами представляют существенный интерес [6–8]. На протяжении десятилетий одним из основных вопросов остается вопрос о влиянии начальных (случайных или детерминированных) возмущений на переход к развитому перемешиванию.

Ориентируясь на задачи ИТС, значительный интерес представляют многослойные системы, в которых возможно не только развитие НРМ, НКГ, но также взаимодействие вторичных волн с зонами перемешивания, развитие вторичных неустойчивостей и т.д. Такой комплексной задачей являются, например, трехслойные газовые системы в ударной трубе. Контактные границы между газами могут иметь различную форму, газы могут иметь различное соотношение плотностей (воздух-SF₆, воздух-He, воздух-Xe), конец трубы может быть закрыт или открыт (прошедшая ударная волна отражается от стенки и повторно взаимодействует с зонами перемешивания). Все эти вариации приводят к существенно различающимся волновым картинам и реализуют различные режимы роста начальных возмущений. Исследованию этих задач посвящено большое число работ [9–19], где экспериментально и численно изучаются различные аспекты реализующегося течения.

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

2. Экспериментальные результаты

Эксперименты проводились на ударной трубе в РФЯЦ-ВНИИЭФ (рис. 1). Техника эксперимента и диагностические методики подробно описаны в предыдущих сходных работах [9, 15, 16]. Ударная труба имеет сечение внутреннего канала 40×40 мм, состоит из камеры высокого давления (драйвера), камеры низкого давления и измерительной секции, состоящей из двух оптически прозрачных отсеков, между которыми в ряде опытов устанавливалась 3D- или 2D-сетка из медной проволочки или из полихлорвиниловых нитей. Контактная граница, на которую накладывалась сетка (КГ₂), имела угол излома 152°. На сетку накладывалась тонкая (34 мкм) полимерная пленка. Размер ячейки 3D-сетки составлял 5×5 мм, толщина медных проволочек — 0,1 мм, толщина нитей ПВХ — 0,06 мм. В 2D-сетке поперечный ряд проволочек отсутствовал, расстояние между проволочками составляло 5 мм.

Рис. 1. Схема постановки опытов.

а) схема ударной трубы; б) фотография измерительной секции с 3D-сеткой 5×5 мм; Д1, Д2 – отметчики времени; КГ₁, КГ₂ – контактные границы

Верхний отсек измерительной секции (КГ₁) отделялся от камеры низкого давления трубы тонкой (3–4 мкм) плоской полимерной пленкой. Камеры высокого и низкого давления трубы разделялись между собой диафрагмой из лавсана толщиной σ $\approx$ 0,1 мм. Камера высокого давления заполнялась сжатым гелием или сжатым воздухом до определенного избыточного давления. В камере низкого давления и в нижнем отсеке измерительной секции находился воздух при атмосферном давлении. Верхний отсек измерительной секции заполнялся либо гелием, либо SF₆ при атмосферных условиях. Нижний торец измерительной секции в ряде опытов был закрыт жесткой горизонтальной стенкой из оргстекла толщиной 20 мм. В вертикальных стенках ударной трубы располагались отметчики времени (пьезокерамические датчики давления Д1, Д2) для определения скорости падающей ударной волны.

Ударная труба работала следующим образом. Саморазрыв диафрагмы, отделяющей камеру высокого от камеры низкого давления, происходил, когда при заполнении драйвера сжатым газом давление превышало критическое. По камере низкого давления (по воздуху) распространялась ударная волна, которая при выходе на КГ₁ инициировала развитие НРМ. Со временем ударная волна достигала КГ₂ и инициировала совместное развитие НРМ и НКГ (из-за излома КГ₂), приводящих к турбулентному перемешиванию газов. Регистрация течения проводилась скоростной видеокамерой шлирен-методом, которая запускалась от электроконтакта.

На рис. 2 и 3 представлены кинограммы некоторых экспериментов. На кинограммах М – число Маха падающей на КГ₁ ударной волны.

Рис. 2. Кинограмма течения в опыте № 1 со слойкой воздух-гелий-воздух (М = 1,25).

1 — гелий; 2 — воздух; 3 – ударная волна; 4 – КГ₂; 5 — ЗТП₁ на КГ₁; 6 — ЗТП₂ на КГ₂; время отсчитывается от прихода падающей ударной волны на КГ₁

Рис. 3. Кинограмма течения в опытах со слойкой воздух-SF₆-воздух.

а) опыт № 6 с 3D-сеткой, М = 1,2; б) опыт № 25 без сетки, М = 1,4.

1 – SF₆; 2 – воздух; 3 – ударная волна; 4 – КГ₂; 5 – ЗТП₁ на КГ₁; 6 – ЗТП₂ на КГ₂; 7 – отраженная от жесткой стенки волна; время отсчитывается от прихода падающей ударной волны на КГ₁

По кинограммам видно, что при проходе УВ из гелия в воздух возмущение от излома КГ₂ не переворачивается (рис. 2), при проходе УВ из SF₆ в воздух – переворачивается (рис. 3). Кроме этого, при отсутствии сетки на КГ₂ из углового возмущения растет ярко выраженный пузырь (локальное возмущение), при наличии сетки пузырь менее выражен. Это говорит о том, что интенсивно развивающаяся из-за сетки зона ТП частично поглощает локальное возмущение. Ширина зоны ТП при наличии сетки больше, чем без сетки, т.е. увеличение шероховатости КГ приводит к увеличению ширины зоны турбулентного перемешивания (ЗТП).

3. Результаты численного моделирования по программам NUT3D и BIC3D

Численное моделирование опыта № 1 в ИПМ РАН было выполнено с использованием двух программных комплексов NUT3D [20-22] и BIC3D [23–26], реализующих две различные численные методики. Методики включают в себя численное решение уравнений невязкой газодинамики в форме Эйлера.

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

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

Общая длина ударной трубы составляла 1730 мм, что при масштабе амплитуды начальных возмущений порядка 0,1 мм и использовании однородной сетки привело бы к числу сеточных элементов порядка 10⁹ и значительным вычислительным ресурсам. По этой причине для трехмерного моделирования был выбран только участок ударной трубы от нижнего торца до координаты по вертикали чуть выше КГ₁, а именно, z =360 мм. В части трубы z > 260 мм задавались параметры течения, соответствовавшие зафиксированной в экспериментах скорости ударной волны. Кроме того, базовый размер ячеек в поперечном направлении (x, y) составлял Δ₁ = 0,4 мм, а в вертикальном направлении (z) – Δ₂ = 0,2 мм. В вертикальной области от 120 мм до 240 мм сетка имела уплотнение с размером ячеек Δ₃ = 0,1 мм. Таким образом, максимальное число ячеек составляло 2,4·10⁷, что позволило существенно ослабить требования к вычислительным ресурсам и выполнить большее число расчетов с различными вариантами начальных возмущений.

Параметры газов соответствовали давлению 1 атм и температуре 20 °C. Скорость ударной волны, инициированной драйвером, составляла 429 м/c.

При прохождении падающей ударной волны через пленку, разделяющую изначально газы, происходит ее разрушение на фрагменты различной величины, которые отвечают за мелкомасштабные начальные возмущения. Разрушение пленки определяется ее внутренним состоянием и в гидродинамическом расчете не моделируется, поэтому возникает необходимость задания реалистичных мелкомасштабных возмущений, имитирующих разрушение пленки. Как известно, от вида начальных возмущений достаточно сильно зависит переход к перемешиванию [27], поэтому желательно задать их наиболее приближенными к эксперименту. Из опытов приблизительно известны диапазон длин волн и амплитуда, определяемые по размерам фрагментов пленки. Следуя работам [13, 28], для имитации разрушения пленки будем задавать случайное возмущение контактных поверхностей в виде:

$\zeta \left(x,y\right)=S{\displaystyle \sum }_{m, n=1}^N\left(\begin{array}{c}{a}_{mn}\cos \left(k_{xm}x\right)\cos \left(k_{yn}y\right)+\\ +b_{mn}\cos \left(k_{xm}x\right)\sin \left(k_{yn}y\right)+\\ +c_{mn}\sin \left(k_{xm}x\right)\cos \left(k_{yn}y\right)+\\ +d_{mn}\sin \left(k_{xm}x\right)\sin \left(k_{yn}y\right)\end{array}\right),$ (1)

где Фурье-коэффициенты a_mn,…,d_mn, выбираются случайным образом на основе нормального распределения, а затем нормируются для получения желаемого среднеквадратичного отклонения амплитуды, S — нормировочный множитель. Волновые векторы определяются выражениями k_xm=2πm/L_x и k_yn= =2πn/L_y, L_x и L_y — поперечные размеры ударной трубы. Суммирование в (1) проводится по допустимому диапазону длин волн, которым будут соответствовать определенные пары волновых чисел m и n. Посредством варьирования спектра начальных возмущений возможно добиться лучшего согласия между численными и экспериментальными данными. Численные расчеты для выбранного опыта были выполнены как в присутствии начальных возмущений на контактных границах, так и без них.

3.2. Физико-математическая модель

В качестве основной физико-математической модели выступала модель многокомпонентной невязкой нетеплопроводной газовой динамики. Соответствующая система уравнений выглядит следующим образом:

$\left\{\begin{array}{c}\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_i}\left(\rho v_i\right)=0\\ \frac{\partial }{\partial t}\left(\rho C_{\alpha }\right)+\frac{\partial }{\partial x_i}\left(\rho C_{\alpha }{v}_i\right)=0\\ \frac{\partial }{\partial t}\left(\rho v_i\right)+\frac{\partial }{\partial x_i}\left(\rho v_i{v}_j+p{\delta }_{ij}\right)=0\\ \frac{\partial }{\partial t}\left(\rho E\right)+\frac{\partial }{\partial x_i}\left(\rho E{v}_i+p{v}_i\right)=\mathrm{0,}\end{array}\right.$ (2)

где ρ — массовая плотность, v_i – компоненты скорости, C_α — массовая концентрация, p – давление, E — удельная полная энергия: E=ε+v_kv_k/2. По повторяющимся индексам подразумевается суммирование. Замыкающие уравнения состояния записываются в модели идеального газа: p_α=(γ_α–1)ρ_αε_α, где ρ_a – фазовые или парциальные плотности веществ, γ_a – показатели адиабаты, ε_a – удельная внутренняя энергия компоненты a. В предположении термодинамического равновесия различных компонент для идеальных газов можно аналитически получить выражение для связи термодинамических параметров смеси, а именно, p=(γ–1)ρε, где

$\gamma =1+\frac{1}{C_V}{\displaystyle \sum }_{\alpha }{C}_{\alpha }{C}_{V\alpha }\left({\gamma }_{\alpha }-1\right)$, $C_V={\displaystyle \sum }_{\alpha }{C}_{\alpha }{C}_{V\alpha }$,

C_V – теплоемкость вещества при постоянном объеме. Использование данной модели физически оправдано и не требует каких-либо дополнительных приближений, так как в качестве исследуемых веществ используются легкие газы. Систему (2) удобно записать в векторном виде, удобном для дальнейшего использования при описании используемых численных схем:

${\partial }_{t}U+{\partial }_i{F}^i\left(U\right)=\mathrm{0,}$ (3)

где U — вектор консервативных переменных, Fⁱ(U) — соответствующие потоковые функции вдоль координатных направлений.

3.3. Численный метод

NUT3D – это трехмерный параллельный программный комплекс, ориентированный на решение задач физики высоких плотностей энергии. Помимо модуля расчета газодинамики (2), NUT3D содержит модули расчета переноса тепла в диффузионном приближении, а также работает с любыми полными термодинамически согласованными уравнениями состояния. Моделирование может вестись в одной из ортогональных систем координат: декартовой, цилиндрической или сферической.

Для получения дискретных аналогов исходных дифференциальных уравнений (2) используется метод конечных разностей. Для этого область моделирования представляется в виде численной ортогональной сетки. Таким образом, дискретная система уравнений может быть записана в следующем виде:

$\begin{array}{c}\frac{({\text{U}}_{\text{ijk}}^{\text{n}+1}-{\text{U}}_{\text{ijk}}^{\text{n}}){\text{V}}_{\text{ijk}}}{\Delta {\text{t}}_{\text{n}}}+\\ +\left({\text{F}}_{\text{i}+1/2\text{jk}}^1{\text{S}}_{\text{i}+1/2\text{jk}}^1-{\text{F}}_{\text{i}-1/2\text{jk}}^1{\text{S}}_{\text{i}-1/2\text{jk}}^1\right)+\dots +\\ +\left({\text{F}}_{\text{ijk}+\frac{1}{2}}^3{\text{S}}_{\text{ijk}+\frac{1}{2}}^3-{\text{F}}_{\text{ijk}-\frac{1}{2}}^3{\text{S}}_{\text{ijk}-\frac{1}{2}}^3\right)={\text{R}}_{\text{ijk}}{\text{V}}_{\text{ijk}},\end{array}$ (4)

где U – вектор консервативных переменных, F – потоки через грани ячеек, S – соответствующие площади граней ячейки, V – объем ячейки. Вычисление потоков через грани ячеек происходит как значений функций, зависящих от решения задачи о распаде произвольного разрыва. В соответствии с работой [29] для определения левых и правых значений вектора переменных на грани используется линейная интерполяция с применением ограничителей антидиффузионных потоков [30]. Решение задачи о распаде разрыва может происходить как точно, так и приближенно. Точное решение задачи о распаде разрыва возможно при использовании УРС идеального газа, а также задействует итерационные методы нахождения корня нелинейного уравнения для давления, вырабатывающегося в результате распада разрыва. По этой причине в расчетах с другими УРС, а также достаточно ресурсоемких, предпочтительнее использовать приближенные методы решения задачи о распаде разрыва. В NUT3D основным методом решения задачи о распаде разрыва является HLLC [31].

В соответствии с методом расщепления по физическим процессам на следующем этапе при необходимости численно решается уравнение теплопроводности:

$\frac{\partial \rho \epsilon \left(\text{T}\right)}{\partial \text{t}}+\nabla {\text{q}}_{\text{T}}=\mathrm{0,}$ (5)

где ${\text{q}}_{\text{T}}=-\kappa \left(\text{T}\right)\nabla \text{T}$. Численное решение уравнения (5) происходит в соответствии с работой [32], в которой предложен итерационный процесс по нелинейности в коэффициенте теплопроводности. Соответствующая численная схема выглядит следующим образом:

$\begin{array}{c}\frac{\left(\epsilon ({\text{T}}_{\text{ijk}}^{\text{s}})+\frac{\partial \epsilon }{\partial \text{T}}({\text{T}}_{\text{ijk}}^{\text{s}})({\text{T}}_{\text{ijk}}^{\text{s}+1}-{\text{T}}_{\text{ijk}}^{\text{s}})-\epsilon ({\text{T}}_{\text{ijk}}^{\text{n}})\right){\rho }_{\text{ijk}}\text{}{\text{V}}_{\text{ijk}}}{\Delta {\text{t}}_{\text{n}}}+\\ +\left({\text{q}}_{\text{i}+1/2\text{jk}}^1{\text{S}}_{\text{i}+1/2\text{jk}}^1-{\text{q}}_{\text{i}-1/2\text{jk}}^1{\text{S}}_{\text{i}-1/2\text{jk}}^1\right)+\dots +\\ +\left({\text{q}}_{\text{ijk}+\frac{1}{2}}^3{\text{S}}_{\text{ijk}+\frac{1}{2}}^3-{\text{q}}_{\text{ijk}-\frac{1}{2}}^3{\text{S}}_{\text{ijk}-\frac{1}{2}}^3\right)=\mathrm{0,}\end{array}$ (6)

где ${\text{q}}_{\text{i}+1/2\text{jk}}^1=-\kappa \left({\text{T}}^{\text{*}}\right)\frac{{\text{T}}_{\text{i}+1\text{jk}}^{\text{s}+1}-{\text{T}}_{\text{ijk}}^{\text{s}+1}}{\Delta {\text{x}}_{\text{ijk}}^1}$, Т^* на грани в соответствии с работой [33] дается выражением

$\text{T}\text{*}=\frac{\kappa \left({\text{T}}_{\text{i}+1\text{jk}}^{\text{s}}\right)\Delta {\text{x}}_{\text{i}-1\text{jk}}^1{\text{T}}_{\text{i}+1\text{jk}}^{\text{s}}+\kappa \left({\text{T}}_{\text{i}-1\text{jk}}^{\text{s}}\right)\Delta {\text{x}}_{\text{i}+1\text{jk}}^1{\text{T}}_{\text{i}-1\text{jk}}^{\text{s}}}{\kappa \left({\text{T}}_{\text{i}+1\text{jk}}^{\text{s}}\right)\Delta {\text{x}}_{\text{i}-1\text{jk}}^1+\kappa \left({\text{T}}_{\text{i}-1\text{jk}}^{\text{s}}\right)\Delta {\text{x}}_{\text{i}+1\text{jk}}^1}$.

Остальные потоки тепла аппроксимируются аналогичным образом. Таким образом, получается система линейных алгебраических уравнений (СЛАУ) относительно вектора переменных Т^s+1.

При расщеплении по координатным направлениям возможно решение исходной СЛАУ (6) последовательным обращением соответствующего количества трехдиагональных матриц методом прогонки. Данный подход может способствовать выделению какого-либо приоритетного направления, поэтому необходимо чередовать последовательность обработки координатных направлений. Также возможно решение исходной СЛАУ (6) каким-либо итерационным методом, например методом сопряженных градиентов с предобуславливателем.

Для вычисления термодинамических параметров в смешанных ячейках используется pT-приближение, предполагающее, что за время гидродинамического шага давление и температура различных компонент смеси успевают выровняться и достигнуть термодинамического равновесия. Также предполагается, что движение всех компонент описывается одной средней скоростью смеси. Кроме того, предполагается, что между компонентами в смешанной ячейке имеется контактная граница. В этом случае массовые и объемные концентрации вводятся выражениями: С_α = m_α /m и f_α = V_α/V. Плотности компонент по определению ρ_α= m_α/V_α = С_α/f_αρ. Таким образом, возникает следующая система уравнений для нахождения равновесного давления, температуры и плотностей компонент (или объемных концентраций) по известным внутренней энергии и плотности смеси после решения задачи о распаде разрыва:

$\left\{\begin{array}{l}\epsilon =\sum C_{\alpha }{\epsilon }_{\alpha }\left({\rho }_{\alpha },T\right);\\ \frac{1}{\rho }=\sum \frac{C_{\alpha }}{{\rho }_{\alpha }};\\ p=p_{\alpha }\left({\rho }_{\alpha },T\right).\end{array}\right.$ (7)

Данная система содержит N_α + 2 уравнений, где N_α — число компонент в ячейке, при необходимости число уравнений может быть уменьшено до N_α.

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

Отметим, что если в качестве независимых термодинамических переменных выбрать р и Т вместо ρ и Т, тогда система, аналогичная (7), содержала бы всего 2 уравнения, что является более приемлемым вариантом с вычислительной точки зрения.

Возможно использование другого подхода к описанию состояния в смешанных ячейках, при котором каждая компонента смеси занимает весь объем ячейки (кинетический подход). В этом случае система уравнений для определения равновесного давления и температуры выглядит следующим образом:

$\left\{\begin{array}{l}\epsilon =\sum {\text{C}}_{\alpha }{\epsilon }_{\alpha }\left({\rho }_{\alpha },\text{T}\right);\\ {\rho }_{\alpha }={\text{C}}_{\alpha }\rho ;\\ p=\sum p_{\alpha }\left({\rho }_{\alpha },\text{T}\right).\end{array}\right.$ (8)

При известных ρ, ε и С_α данная система требует решения только уравнения для нахождения равновесной температуры.

BIC3D – это параллельный программный комплекс для численного решения трехмерных эволюционных гиперболических систем законов сохранения по бикомпактным схемам из [23–26].

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

$\Omega ={\Omega }_1^x\times {\Omega }_1^y\times {\Omega }_1^z=\{(x_j,y_k,z_l)\};$

${\Omega }_1^d=\left\{d_0<d_{1/2}<d_1<d_{3/2}<d_2<\dots <d_{N_d}\right\};$

 $h_{d,i}=d_{i+1}-d_i$ $–$ шаг по d, $i=\mathrm{0,} \mathrm{1,}\dots , N_d-1.$ (9)

На временной оси вводятся слои t_n, n =0, 1, ...; t =τ_n= t_n+1 – t_n — шаг по t. Численное решение представляется в виде сеточной функции $U_{i,j,k}^n$, значения которой аппроксимируют значения точного решения U(x,y,z,t) в узлах сетки Ω в моменты времени t =t_n.

В бикомпактных схемах из [23, 24] используется локально-одномерное расщепление по пространству. Помимо него, осуществляется потоковое расщепление Лакса–Фридрихса–Русанова (ЛФР). Совместное применение этих расщеплений означает, что при каждом переходе со слоя на слой трехмерная система (2) разбивается на шесть одномерных подсистем вида

${\partial }_{t}U+{\partial }_d{F}^{d,\pm }\left(U\right)=\mathrm{0,}$ (10)

где $F^{d,\pm }(U)=0.5F^d(U)\pm C_2^{d}U$ (формула для $C_2^d$ будет приведена позднее). Одномерные подсистемы (2) решаются при помощи одномерных бикомпактных схем вдоль линий сетки – прямых, проходящих через узлы сетки Ω параллельно координатным осям. Обозначим оператор послойного перехода всякой одномерной бикомпактной схемы для (2) как $S_B^{d,\pm }(\tau ).$ Введем операторы $S_B^d(\tau )=S_B^{d,-}(\tau )S_B^{d,+}(\tau ).$ Результирующая локально-одномерная бикомпактная схема для исходной системы (2) имеет оператор послойного перехода

$S_B\left(\tau \right)=S_B^x\left(\frac{\tau }{2}\right)S_B^y\left(\frac{\tau }{2}\right)S_B^z\left(\tau \right)S_B^y\left(\frac{\tau }{2}\right)S_B^x\left(\frac{\tau }{2}\right).$ (11)

Заметим, что точность пространственного расщепления повышена по методу Марчука–Стрэнга.

Аппроксимация по времени в бикомпактных схемах строится при помощи чисто неявных и неявно-явных методов Рунге–Кутты (РК). В BIC3D реализованы методы РК обоих типов от первого до четвертого порядков включительно. Временные дискретизации четных порядков лучше подходят для счета гладких решений, нечетных – для счета разрывных решений. Неявно-явные бикомпактные схемы экономичнее, так как в них исключены ньютоновские итерации; чисто неявные бикомпактные схемы точнее и генерируют меньшие осцилляции у фронтов сильных разрывов.

В расчетах, результаты которых представлены в следующем разделе, операторы $S_B^{d,\pm }(\tau )$ отвечали неявно-явной схеме BiC4-IMEX3 из [24]. Эта схема имеет четвертый и третий порядки аппроксимации по пространству и времени соответственно.

Для подавления эффекта Гиббса в бикомпактных схемах повышенного порядка аппроксимации по времени применяется метод консервативной монотонизации [25]. Он имеет один настраиваемый параметр C₁≥ 0. При C₁= 0 монотонизация отсутствует, при C₁→∞ она максимальна (от численного решения в каждой ячейке остается лишь интегральное среднее). Метод [25] основан на сравнении решения высокоточной схемы (B, в данном случае бикомпактной) с решением монотонной схемы-партнера (A). В BIC3D в качестве схемы A выбран и реализован локально-одномерный “явный уголок” с расщеплением ЛФР. Оператор послойного перехода этой схемы имеет вид:

$S_A\left(\tau \right)=S_A^z\left(\tau \right)S_A^y\left(\tau \right)S_A^x\left(\tau \right).$ (12)

Здесь $S_A^d\left(\tau \right)=S_A^{d,-}\left(\tau \right)S_A^{d,+}\left(\tau \right)$, где $S_A^{d,\pm }\left(\tau \right)$ – операторы послойного перехода одномерного “явного уголка” для (2).

Коэффициенты $C_2^d$ и шаг $\tau$ задаются перед каждым переходом со слоя на слой по формулам

$\begin{array}{c}{C}_2^d=\left(0.5+\delta \right)V_{\text{max}}^d;\\ \tau =\frac{\kappa }{1+\delta }\underset{d}{\min }\frac{{\text{min}}_i{h}_{d,i}}{V_{\text{max}}^d}; \\ V_{\text{max}}^d=\underset{\Omega }{\max }{V}^d(U_{j,k,l}^n),\end{array}$ (13)

где V^d(U) – мажоранты для модулей собственных значений матриц $A^d(U)={\partial }_U{F}^d(U)$, δ > 0 – коэффициент запаса, κ > 0 – число Куранта.

3.4. Программная реализация

BIC3D реализован на языке C++ в рамках технологии MPI и предназначен для использования на многопроцессорных вычислительных системах. В соответствии со структурой численного подхода (см. предыдущий раздел), а именно, использование локально одномерного расщепления по координатным направлениям, в качестве идеологии распараллеливания используется геометрический параллелизм. В этом случае в пределах одной стадии такого расщепления вычисления на параллельных линиях сетки не связаны между собой, и потому могут осуществляться одновременно. Метод консервативной монотонзации из [25] также легко распараллеливается: для построения монотонизированного решения в каждом узле нужны данные лишь от ячеек, которым он принадлежит.

Пусть число процессов равно n_p. Разобьем всю совокупность ячеек сетки на n_p горизонтальных слоев равной толщины M_z=N_z/n_p (в ячейках) и n_p вертикальных слоев равной толщины M_x=N_x/n_p. Ради простоты записи мы полагаем, что N_z, N_x делятся нацело на n_p. Процесс ранга $r=\overline{\mathrm{0,}{n}_p-1}$ имеет дело не со всей сеткой Ω, а лишь с ее подмножествами

${\omega }_r={\Omega }_1^x\times {\Omega }_1^y\times \left\{z_l:0\le l-r{M}_z\le M_z\right\};$

${{\omega }^{\prime }}_r=\left\{x_j:0\le j-r{M}_x\le M_x\right\}\times {\Omega }_1^y\times {\Omega }_1^z.$ (14)

Соответственно, процесс ранга r оперирует двумя сеточными функциями (массивами) решения, определенными на сетках ${\omega }_r$, ${{\omega }^{\prime }}_r$.

Опишем алгоритм перехода со слоя t_n на слой t_n+1 с точки зрения процесса ранга r.

1.       Найти $\underset{{\omega }_r}{\max }{V}^d(U_{j,k,l}^n)$. Если r ≠0, то отправить эти максимумы процессу ранга 0 и в ответ получить величины $V_{\text{max}}^d$. Если r =0, то собрать данные от остальных процессов, вычислить величины $V_{\text{max}}^d$ и разослать их остальным процессам. Вычислить $C_2^d$ и τ по формулам (4).

2.       Рассчитать на ω_r промежуточное решение

$Y^B=S_B^y\left(\frac{\tau }{2}\right)S_B^x\left(\frac{\tau }{2}\right)U^n.$ 

Действие, например, $S_B^{x,+}(\tau /2)$ на $U^n$ реализуется так: с помощью одномерной бикомпактной схемы для всех $0\le k\le N_y$, $0\le l-r{M}_z\le M_z$ решить подсистему (2) с индексами (x, +) вдоль прямой y = y_k, z = z_l при половинном временном шаге, используя в качестве начальных данных Uⁿ.

3.       Разбить массив решения Y^B на ω^r сообразно разбиению множества ячеек на вертикальные слои; полученные части этого массива разослать по соответствующим процессам (рис. 1); собрать данные, присланные другими процессами. Как результат, у данного процесса в распоряжении окажется все то же решение Y^B, но теперь уже на ${{\omega }^{\prime }}_r$.

4.       Рассчитать на ${{\omega }^{\prime }}_r$ промежуточное решение ${\tilde{Y}}^B=S_B^z(\tau )Y^B.$

5.       Выполнить операцию, обратную (симметричную по x, z) операции из шага 3. У данного процесса получится решение ${\tilde{Y}}^B$ на ω^r.

6.       Рассчитать на ω^r финальное решение схемы B на слое t_n+1: $U^B=S_B^x\left(\frac{\tau }{2}\right)S_B^y\left(\frac{\tau }{2}\right){\tilde{Y}}^B.$

7.       Получить на ω^r вспомогательное решение U^A, повторяя шаги 2–5 для схемы A, учитывая при этом формулу (3) для ее оператора послойного перехода.

8.       Провести монотонизацию решения U^B при помощи решения U^A по методу из [22]. При исполнении этой операции необходимо обменяться с процессами r ±1 данными с плоскостей l=rM_z, l= =(r+1)M_z.

9.       Искомое решение $U^{n+1}$ на ω^r найдено, переход на слой t_n+1 завершен.

Барьерные синхронизации процессов выполняются в начале каждого перехода со слоя на слой и перед операциями обменов. Ясно, что при C₁ = 0 из алгоритма исключаются шаги 7, 8.

В работе [26] было проведено исследование BIC3D на эффективность распараллеливания: она составляет от 70% при условии $n_p\le ({\min }_d{N}_d)/2.$

В завершение дадим практические рекомендации по работе с BIC3D. При использовании комплекса пользователь задает:

• Физическую модель – функции F^d(U) и V^d(U). Решение задачи Римана (Riemann solver) специфицировать не нужно из-за глобального расщепления ЛФР.
• Пространственную сетку Ω через задание сеток ${\Omega }_1^d$ (как массивов узлов).
• Начальные условия в виде функции для вычисления начальных значений либо уже готового файла данных, содержащего эти значения (очевидно, размеры файла должны соответствовать сетке Ω).
• Типы условий на всех границах расчетной области. Поддерживаются: условия первого рода (для них необходимо дополнительно задать функции для вычисления граничных значений), условия отражения, условия свободного выхода и периодические условия.
• Параметры δ, κ, C₁. Рекомендуемые значения: C₁=0 и любое κ при счете гладких решений, C₁=10, κ≥0.6 при счете разрывных решений (диссипация в бикомпактных схемах нечетного порядка растет с увеличением κ), δ=0.5. При κ>1 шаг по времени в схеме A автоматически дробится так, чтобы она сохраняла устойчивость.

   

  

  Рис. 4. Рассылка решения с горизонтального слоя (n_p=4, r =2).

   

  Численный код NUT3D в последовательном варианте реализован на языке C. Параллельная версия, использовавшаяся для моделирования в данной работе, была получена с помощью набора ин­струментальных средств, созданных в ИПМ им. М.В. Кел­дыша РАН и упрощающих процесс распараллеливания. Применение этих инструментов состоит из нескольких этапов. На первом этапе программа на языке Fortran или C подается на вход системе автоматизированного распараллеливания (далее системе SAPFOR [34, 35]). Блок анализа системы переводит программу во внутреннее представление. Программа анализируется с помощью различных алгоритмов статического анализа, в результате которого выявляются особенности программы, влияющие на ее распараллеливание. Затем на втором этапе происходит подбор подходящих схем распараллеливания программы. Отбор подходящих схем распараллеливания осуществляется на основе прогнозируемых характеристик эффективности выполнения параллельной программы. Система может предлагать пользователю на выбор разные схемы распараллеливания, сопровождая их прогнозируемыми характеристиками эффективности и причинами отказа от распараллеливания фрагментов программы. Пользователь определяет значимые фрагменты программы и исследует причины отказа от их распараллеливания. Возможно уточнить результаты автоматического анализа (например, используя дополнительные спецификации) и повторить второй шаг этап, и/или преобразовать текст последовательной программы и начать с первого. При успешном выполнении этого этапа система по отобранным пользователем схемам распараллеливания и внутреннему представлению программы генерирует текст параллельной программы в модели DVMH на языке Fortran-DVMH или C-DVMH.

  Модель программирования DVMH [36] позволяет разрабатывать параллельные программы для кластеров, в узлах которых помимо универсальных многоядерных процессоров установлены ускорители и сопроцессоры.

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

  Для отображения DVMH-программы на параллельную ЭВМ используются специальные компиляторы, которые входят в состав DVM-системы.

  DVMH-компилятор преобразует исходную программу в параллельную программу на языке Fortran или C с вызовами функций системы поддержки параллельного выполнения DVMH-программ (библиотеки Lib-DVMH). Система поддержки реализована на основе стандартных технологий параллельного программирования MPI (для взаимодействия между узлами кластера), OpenMP (для распределения вычислений между ядрами узла), CUDA (для распределения вычислений и отображения данных на графические ускорители).

  3.5. Результаты расчетов и сравнение с экспериментом

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

   

  

  Рис. 5 (начало).

   

  

  Рис. 5 (продолжение).

   

  

  Рис. 5 (продолжение).

   

  

  Рис. 5 (продолжение).

   

  

  Рис. 5 (окончание). Экспериментальные фотографии (а) и численные шлире

   

  В приведенной на рис. 5 картине течения (t =0,052 мс) на экспериментальной фотографии в кадре нет верхней (плоской) контактной поверхности. При этом видна ударная волна, прошедшая в гелий, положение фронта которой совпадет и в эксперименте, и в численных расчетах. В следующий приведенный момент времени (t =0,132 мс) падающая ударная волна уже прошла вторую контактную границу и вышла из гелия, изменив форму фронта. Здесь уже наблюдаются различия в координате фронта этой ударной волны – в численных расчетах она распространяется быстрее. Далее (t =0,232 мс) данная тенденция сохраняется. На численных шлирен-изображениях наблюдается ударная волна, отраженная от КГ₂ и прошедшая через КГ₁ (z =224 мм). В t =0,352 мс в кадре появляется ­возмущенная КГ₁, однако расчетное ее положение существенно отличается от экспериментального, отставая от него. Также отставание начинает наблюдаться и для КГ₂, движение которой в более ранние моменты времени совпадало с экспериментальным. Финальный момент времени t =0,532 мс лишь подтверждает наметившуюся тенденцию. На основании рис. 5 можно отметить следующие обстоятельства. Во-первых, расчеты, выполненные с использованием двух различных методик (NUT3D и BIC3D), хорошо согласуются между собой. Во-вторых, в эксперименте скорость движения контактных границ выше, чем в расчетах. Это может быть связано с неточными данными по положению пленки, разделяющей воздух и гелий в измерительной секции трубы, и скорости падающей ударной волны. Кроме того, в расчетах не учитываются особенности взаимодействия падающей ударной волны с пленкой, которые могут приводить к ее выгибанию и смещению перед разрушением.

  В следующей серии расчетов на обеих контактных границах были заданы случайные начальные возмущения. Длина волны возмущений задавалась в интервале от λ_min=1 до λ_max=3 мм, среднеквадратичное отклонение составляло 0,1λ_min. На рис. 6 приведены шлирен-изображения, полученные в расчетах с использованием NUT3D и BIC3D.

  На основании рис. 6 можно отметить следующее. Наличие возмущений на контактных границах незначительно влияет на скорость их движения – положения в различные моменты времени идентичны расчетам без возмущений (рис. 5). Характер роста возмущений на КГ₁ и КГ₂ отличается, во-первых, в связи с дополнительным развитием НКГ на КГ₂, а во-вторых, вследствие инверсии фаз начальных возмущений на КГ₁ при прохождении исходной падающей ударной волны. Сформировавшиеся в результате развития начальных возмущений зоны перемешивания не являются идентичными. На КГ₂ можно отчетливо рассмотреть вихреобразные структуры, а именно дорожку вихрей, характерную для развития неустойчивости Кельвина-Гельмгольца. Также можно отметить, что развитие возмущений на этой контактной границе происходит медленнее, чем на КГ₁.

   

  

  Рис. 6. Численные шлирен-изображения в различные моменты времени: (а) NUT3D, (б) BIC3D.

   

  4. Верификация численных кодов ЭГАК и k-ε модели

  В данном разделе представлены результаты проверки численных кодов ЭГАК и k-ε модели на эксперименте № 6.

  Методика ЭГАК является лагранжево-эйлеровой, в которой аппроксимация уравнений многокомпонентной газодинамики производится в два этапа на неподвижной прямоугольной, как правило, квадратной счетной сетке [37].

  На лагранжевом этапе вычислений используется схема типа “крест”. Все компоненты среды выделяются своими термодинамическими параметрами и объемной долей. Для замыкания уравнений газодинамики в смешанных ячейках, содержащих два и более компонента, используются несколько методов, основанные на разных предположениях относительно термодинамического состояния компонентов [38].

  На эйлеровом этапе в чистых (однокомпонентных) ячейках для аппроксимации уравнения адвекции используется метод PPM [39] третьего порядка аппроксимации. Для определения потоков из смешанных ячеек используется метод концентраций (VоF) [40].

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

  Для моделирования турбулентного перемешивания в методике используется как k-ε модель, так и прямое моделирование на подробной сетке без использования каких-либо моделей турбулентности. В настоящей работе использовался подход, основанный на решении уравнений невязкой газодинамики в форме Эйлера.

  4.1. Математическая постановка задачи

  Начальная геометрия эксперимента взята из рис. 1. На КГ₂ в эксперименте была задана двухмерная или трехмерная сетка из медной или лавсановой проволоки с λ=0,5 см. В расчетах эта сетка была неподвижной и моделировалась с использованием несжимаемого компонента. При этом в расчетах диаметр проволоки задавался равным d_расчет=d +2s, где d — реальный диаметр проволоки, s — толщина натянутой на сетку пленки. Это сделано исходя из предположения, что пленка после разрыва в процессе течения прилипает к проволоке и ее фактическая толщина увеличивается. В двух расчетах сама пленка не задавалась, в третьем расчете она задавалась в виде разорванных кусков толщиной 0,005 см и шириной 0,03 см, расстояние между кусочками составляло 0,01 см.

  В расчетах используется следующая система единиц: г, см, мс = 10⁻³ с, Т, К, давление в атм. На КГ₁ задавалось условие свободного протекания. На боковых стенках задавалось условие “жесткая стенка”. КГ₁ была возмущена в виде синусоиды с длиной волны λ=0,03 см и с амплитудой возмущений а = 0,01 см.

  Уравнение состояния газов:

  $\epsilon =c_V\cdot T=\frac{p}{(\gamma -1)\rho }$,

  здесь γ — постоянная адиабаты. В расчетах с пленкой для пленки использовался УРС Ми—Грюнайзена с ρ₀=0,92, c₀=2,65, n =4,5, γ=2,35, c_V =1, p₀=1.

  Начальные данные задачи зависели от использованных в опытах газов. Отметим, что в области “газ 1” (воздух) задавался соответствующий газ с параметрами за фронтом ударной волны, которые определялись следующим образом.

  Начальное давление в области “газ 1”

  $p_1=\frac{2\cdot {\rho }_0\cdot u_{\text{óâ}}^2-({\gamma }_1-1)\cdot {ð}_0}{({\gamma }_1+1)}$. (15)

  Здесь скорость ударной волны u_ув берется из эксперимента.

  Далее находим начальную плотность в “газ 1”

  ${\rho }_1={\rho }_0\cdot \frac{({\gamma }_1+1)\cdot {ð}_1+({\gamma }_1-1)\cdot {ð}_0}{({\gamma }_1+1)\cdot {ð}_0+({\gamma }_1-1)\cdot {ð}_1}$. (16)

  И наконец, находим начальную скорость в “газ 1”

  $u_z=-\sqrt{(p_1-p_0)\cdot \left(\frac{1}{{\rho }_0}-\frac{1}{{\rho }_1}\right)}$.   (17)

  Начальные данные опыта № 6: в области “газ 1”: ρ₁₀ = 0,0017, u_z0=21,62, р₁₀=1,5844, γ₁=1,4; в области “газ 2”: ρ_{SF6 0}= 6,5·10⁻³, u₀=0, р₂₀=1, γ₂=1,0945; в области ниже КГ₂ ρ₀=1,25·10⁻³, u₀=0, р₀=1, γ₁=1,4.

  Проведены три расчета на счетной сетке с h =0,005: расчет 1 – без пленки и k-ε модели; расчет 2 – без k-ε модели с кусками пленки, заданными как указано выше; расчет 3 – без пленки с k-ε моделью. В расчете с k-ε моделью на КГ₂ задавались начальные затравочные значения турбулентной энергии и скорости ее диссипации: k =10⁻⁴, ε=10^-5.

  Численное моделирование задачи производилось в 2D-постановке по методике ЭГАК [37].

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

  Ниже приводится сравнение с экспериментами (с двумерной 2D- и с трехмерной 3D-сеткой) указанных вариантов двумерных расчетов. Эксперименты представлены шлирен-фотографиями, расчеты – распределением градиента плотности, по сути, это расчетная шлирен-картина (рис. 7–10).

   

  

  Рис. 7. Картины течения: расчеты t = 491, эксперимент с 3D-сеткой t = 491,6, эксперимент с 2D сеткой t = 481.

   

  

  Рис. 8. Картины течения: расчеты t = 571, эксперимент с 3D-сеткой t = 571,6, эксперимент с 2D-сеткой t = 565,1.

   

  

  Рис. 9. Картины течения: расчеты t = 732, эксперимент с 3D-сеткой t = 731,6, эксперимент с 2D-сеткой t = 732,1.

   

  

  Рис. 10. Картины течения: расчеты t = 812, эксперимент с 3D-сеткой t = 811,6, эксперимент с 2D-сеткой t = 815,6.

   

  Отметим некоторые различия, имеющиеся в экспериментальных данных. Имеется некоторое опережение движения волн в эксперименте с 2D-сеткой. На рис. 7 видно, что ударная волна пришла на сетку на более раннее время по сравнению с 3D-сеткой. Такое небольшое отличие может быть следствием изменения сетки или же это вследствие экспериментальной погрешности в определении скорости ударной волны, движущейся по ударной трубе. Как бы там ни было, в этой ситуации от расчетов правомерно ожидать качественного согласия с экспериментальными данными, не претендуя на большую точность по времени. Кроме того, имеются отличия в интенсивности перемешивания и скорости распространения ЗТП вдоль трубы, что хорошо видно на рис. 10.

  Отметим, что на всех расчетных рисунках приводятся распределения градиентов плотности в диапазоне от 0 до 0,01, соответственно этому цвет меняется от синего до красного. Таким образом, красные оттенки в расчетах дают представление о процессах перемешивания, в частности о зоне турбулентного перемешивания.

  Из сравнения расчетных и экспериментальных распределений можно сделать следующие выводы. В целом результаты расчетов согласуются между собой и с экспериментальными данными, при этом расчеты лучше согласуются с экспериментом с 2D-сеткой, что не удивительно, так как они моделируют именно этот опыт (расчеты двумерные). Наличие проволочной сетки в расчете 1 чувствуется до конца моделирования, в расчетах хорошо видны струйки с большими градиентами, совпадающими с положением сетки. Менее ярко этот эффект проявляется в расчете 2, наличие фрагментов пленки хаотизирует течение и эти струйки со временем размываются. Еще более выделяется расчет с k-ε моделью, в котором наличие диффузии сглаживает этот эффект с самого начала течения, но он все равно заметен (рис. 10).

  Однако имеются и более существенные отличия расчетов от экспериментов. Во всех расчетах есть небольшое отставание ударной волны, приходящей на место расположения проволочной сетки (рис. 7), по сравнению с экспериментальными данными, на рис. 9 видно, что расчетная ударная волна опережает экспериментальную. Наиболее согласующиеся с экспериментом данные получены в расчете 2 с учетом фрагментов пленки. В частности, на рис. 9 в расчете 2 отраженная от пленки волна (которая распространяется вверх по трубе) заметно ближе к экспериментальному положению. Еще более заметно, что в расчете 2 ЗТП на границе воздух–SF₆ согласуется с экспериментальными данными, в то время как в расчете 1 она практически отсутствует, а имеются лишь завитушки вследствие проявления сдвиговой неустойчивости, а в расчете 3 ЗТП имеется, однако она значительно тоньше. Таким образом, пленка играет важную роль в данном течении и для корректного численного моделирования необходимо ее учитывать.

  Отметим, что для 3D-сетки надо бы провести 3D-моделирование, однако подобный расчет требует неоправданно больших ресурсов компьютера. На качественном уровне результаты вполне адекватны и при 2D-моделировании.

  5. Верификация численной методики МИМОЗА

  Моделирование опыта № 1 выполнялось в постановке, приведенной на рис. 1. В трехслойной газовой системе центральный слой (между контактными границами КГ₁ и КГ₂) заполнялся легким газом (He). Вне центрального слоя располагался воздух. Все газы находились при атмосферном давлении.

  УВ формировалась на левой границе и двигалась направо с числом Маха М=1,3. Граничные условия: на верхней и нижней границах ставились периодические условия, за фронтом УВ газодинамические параметры соответствовали величинам за фронтом скачка.

  Выполнено два расчета. В первом расчете на обеих КГ задавались двухмодовые (а₁, а₂≠0, а₃=0) начальные возмущения, во втором – трехмодовые (а₁, а₂, а₃ ≠ 0) возмущения по соотношению

  $D(y)=a_1\sin (k_{1}y)+a_2\sin (k_{2}y)+a_3\sin (k_{3}y).$

  Амплитуды и длины волн составляли а₁=0,01 см, λ₁=3 см, а₂=0,01 см, λ₂ = 0,4 см, а₃=0,002 см, λ₃=0,01 см.

  5.1. Методика расчета

  Математическое моделирование опыта №1 выполнялось в 2D-постановке по методике МИМОЗА с использованием уравнений Эйлера без привлечения каких-либо моделей учета ТП (ILES моделирование). Расчетная методика основана на лагранжево-эйлеровой стратегии и выделении веществ концентрациями (подробнее см. [41–44]). Такой подход является эффективным при моделировании задач механики сплошной среды с большими деформациями.

  Расчет счетного шага состоит из двух этапов: на первом этапе выполняется интегрирование уравнений Эйлера, записанных в лагранжевых координатах, на втором этапе производится пересчет полученных сеточных значений на первоначальную квадратную сетку [42]. Пересчет величин осуществляется при помощи алгоритма, основанного на расщеплении по координатным направлениям и использовании одномерного алгоритма повышенного порядка точности.

  На лагранжевом этапе расчета границы ячеек сетки перемещаются со скоростью вещества, массы ячеек не изменяются. Интегрирование системы уравнений выполняется на разнесенной разностной сетке. Термодинамические параметры задачи относятся к центру счетной ячейки, координаты и компоненты скорости – к узлам. Используется полностью консервативная разностная схема “предиктор–корректор”, аналогичная [44]. С целью предотвращения размытия КГ между веществами применяется алгоритм их выделения с помощью концентраций. Для подавления паразитических осцилляций численного решения в окрестности больших градиентов газодинамических величин вводится искусственная вязкость, являющаяся суммой квадратичной и линейной вязкостей. В смешанных ячейках (ячейках, содержащих несколько веществ) давление вычисляется покомпонентно, а затем усредняется с учетом объемных концентраций веществ.

  Расчет многокомпонентной сплошной среды осуществляется посредством двух методов отслеживания КГ материалов: метода концентраций – VolumeofFluid (VoF) [43] и метода моментов концентраций – MomentofFluid (MoF) [45]. Особенностью метода VoF является зависимость алгоритма от объемных концентраций материала в соседних счетных ячейках. Другая особенность VoF-методов заключается в том, что при наличии в смешанной ячейке более двух компонентов результат реконструкции КГ зависит от порядка обработки веществ в алгоритме.

  В методе MoF, зная положение центра массы вещества и занимаемую им долю объема в ячейке, путем решения задачи минимизации определяется КГ вещества в виде прямой линии без привлечения данных из соседних счетных ячеек. Метод MoF имеет второй порядок точности, что позволяет точно реконструировать линейные КГ. Метод MoF позволяет точнее рассчитывать потоки объемов компонент на этапе адвекции и позволяет дольше удерживать целостность КГ материалов в процессе счета.

   

  

  Рис. 11. Картины течения: расчеты t = 1000, эксперимент с 3D сеткой t = 1071,6 (для эксперимента с 2D-сеткой данных нет).

   

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

  Результаты моделирования динамики границ и зон перемешивания, полученные с использованием методов VoF и MoF, дали близкие результаты. На рис. 12 сопоставляются результаты экспериментов и расчетов, полученных методом VoF. Рассчитанные по методике МИМОЗА зоны перемешивания соприкасающихся газов на КГ₁ и КГ₂ изображены на рисунке зеленым цветом.

   

  

  Рис. 12. Сравнение расчетной и экспериментальной динамики границ и зон перемешивания, слойка воздух-гелий-воздух, t ≈530 мкс.

   

  Анализ полученных результатов показал, что:

• ширина области перемешивания при трехмодовых начальных возмущениях заметно меньше, чем при двухмодовых (наиболее ярко это наблюдается на КГ₁);
• расчетная скорость КГ₁ несколько меньше, чем в эксперименте.      

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

  Проведены численные расчеты динамики роста малых возмущений на контактных границах газов в трехслойных системах воздух–He–воздух и воздух–SF₆–воздух вследствие развития НРМ и НКГ. Эти данные сопоставлены с данными натурных экспериментов, выполненных в ударной трубе. Получено удовлетворительное согласие расчетных и экспериментальных результатов по динамике зоны турбулентного перемешивания. Некоторые временные отличия могут быть связаны с экспериментальными погрешностями измерения скорости падающей на контактную границу ударной волны.

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

  Расчеты ИПМ им. М.В. Келдыша РАН выполнены с использованием вычислительных ресурсов Суперкомпьютерного центра коллективного пользования.

  П.А. Кучугов выражает благодарность разработчикам системы SAPFOR и DVMH за реализацию параллельной версии программы NUT3D.

  Работа выполнена в рамках научной программы Национального центра физики и математики по Государственному контракту № Н.4ц.241.4Д.23.1085.