Models for description of subsonic flows with premixed turbulent combustion in channels

全文:

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

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

В отличие от турбулентных течений без горения, метод крупных вихрей (LES) не гарантирует улучшение предсказания турбулентного горения по сравнению с расчетами на базе уравнений Рейнольдса (RANS) (см., напр., [1]). Дело в том, что в турбулентном горении существенную роль играет молекулярное смешение топлива с окислителем и молекулярная диффузия тепла, которые происходят на масштабах мельчайших турбулентных вихрей. Поскольку LES не позволяет опуститься до столь мелких масштабов, то даже в LES приходится использовать приближенные (полуэмпирические) модели микросмешения и горения, что в значительной степени понижает превосходство LES-расчета над RANS-расчетами, приводит к падению точности результатов и необходимости настройки эмпирических констант подсеточной модели. Правда, сохраняется важное преимущество LES-расчета перед RANS-расчетом – намного более высокая точность определения локальных условий горения, что связано с использованием гораздо более подробных расчетных сеток. Данная статья ограничивается моделями, предназначенными к применению в рамках подхода RANS.

Для замыкания системы уравнений RANS или LES необходимы три модели: модель турбулентности (в случае LES – модель подсеточной турбулентности), модель взаимодействия турбулентности и горения (Turbulence–Combustion Interaction, TCI) и модель химической кинетики. К сожалению, физико-математические модели, описывающие корректно многомасштабные физические процессы, протекающие в турбулентном пламени (смешение компонентов, диффузию тепла, воспламенение, стабилизацию и срыв горения), до сих пор не созданы. Качество моделирования явления в значительной степени зависит от правильного выбора перечисленных трех моделей.

В настоящей работе основное внимание будет уделяться проблемам моделирования TCI и выбора подходящей модели химической кинетики. В качестве модели турбулентности будут использоваться двухпараметрические дифференциальные модели класса k - ω, использующие специальную функцию для плавного перехода от пристеночной (низкорейнольдсовой) версии модели к модели, рассчитанной на описание свободной турбулентности. Это естественный выбор при описании течений в каналах, где есть области и пристеночной, и свободной турбулентности. Дополнительными полезными свойствами моделей класса  k - ω являются их хорошие свойства в окрестности твердых стенок (модели этого класса имеют хорошую асимптотику у стенок без использования пристеночных коррекций [2]), а также то, что они позволяют оценить локальное значение характерного интегрального временного масштаба турбулентных пульсаций скорости τ = 1 / ω. Отметим, что первая дифференциальная модель этого класса (и вообще первая в истории дифференциальная модель турбулентности) была сформулирована еще А.Н. Колмогоровым [3].

Наибольшие трудности возникают с описанием взаимодействия турбулентности и горения (TCI). Детальное обсуждение проблем турбулентного горения можно найти в книгах [4–8].

Существует два основных канала TCI. 1-й канал TCI связан с влиянием турбулентности на средние скорости химических реакций, 2-й – с влиянием горения на турбулентные потоки тепла и массы компонент реагирующей смеси. Существуют и другие пути взаимодействия турбулентности и горения [9], но они имеют косвенный характер (т.е. представляют собой цепочки взаимодействий различных факторов, которые приводят, в конечном счете, к влиянию турбулентности на горение или обратно).

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

В практических расчетах 1-й канал TCI часто не учитывается, и используется т.н. квазиламинарное приближение, при котором средние скорости реакций вычисляются по обычным формулам, в которые подставляются средние параметры течения. Однако квазиламинарное приближение часто приводит к большим ошибкам, поскольку при турбулентном горении пульсации параметров могут быть сопоставимы со средними величинами и даже превосходить их на порядок величины. Поэтому к настоящему времени предложено много способов описания 1-го канала TCI. Среди этих способов стоит выделить следующие классы:

1. статистические методы (метод моментов, метод функции плотности вероятности – ФПВ) – наиболее формальный подход, учитывающий статистические характеристики турбулентного течения с горением, но не опирающийся на физические представления о структуре турбулентного пламени [10–14];
2. модели микроламинарных пламен (флеймлетов) как для неперемешанного, так и для предварительно перемешанного горения – подход, применимый в случае, когда характерное время тепловыделения мало по сравнению с характерным временем молекулярного смешения на уровне мельчайших турбулентных вихрей [15–18]; существуют расширения этого подхода на другие случаи, связанные с введением переменной прогресса реакции (например, [19]), но они опираются на много сильных допущений и менее надежны;
3. модели реактора частичного перемешивания (PaSR – Partially Stirred Reactor) [20–23]. В этих моделях предполагается, что при высоких числах Рейнольдса горение протекает в основном в т.н. “тонких структурах”, связанных с мелкомасштабной турбулентностью. В моделях класса PaSR тонкие структуры рассматриваются как реакторы, в которых непрерывно протекает горение. В работах [24, 25] представлен успешный опыт применения моделей класса PaSR к описанию утолщенных пламен, в которых толщина фронта пламени в несколько раз превышает микромасштаб Колмогорова;
4. фронтальные модели для предварительно перемешанного горения [26–28], основанные на оценке скорости турбулентного пламени. Стоит отметить, что первые оценки этой скорости для различных режимов турбулентного горения были получены еще К.И. Щелкиным [29] и А.Г. Прудниковым [30]. Во фронтальныx моделях, как правило, делается предположение, что в случае сильной турбулентности, типичном для практических приложений, горение происходит в мгновенных тонких реакционных зонах, структура которых слабо отличается от реакционных зон в ламинарном пламени. Такая картина течения была получена в ряде экспериментальных данных как для предварительно перемешанного горения [31–33], так и непремешанного горения [34], а также средствами прямого численного моделирования турбулентности [35].

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

Менее исследован 2-й канал TCI. Нередко считают, что разработка моделей 2-го канала TCI важна в RANS-моделях течения (на основе уравнений Рейнольдса) и не важна при использовании вихреразрешающих подходов к моделированию турбулентных течений, т.к. основную роль во 2-м канале TCI играет крупномасштабная турбулентность, которая в вихреразрешающих подходах описывается напрямую. Однако при описании течений в каналах и обтекания тел приходится применять гибридные вихреразрешающие методы [36], которые в пристеночной области течения переходят в моделирование на базе нестационарных уравнений RANS. В этом случае в пристеночной области учет 2-го канала TCI оказывается существенным.

Чаще всего используется гипотеза об изотропном характере переноса скалярных величин (энтальпии и массовых долей компонент смеси), и тогда 2-й канал TCI учитывают путем введения переменных турбулентных чисел Прандтля и Шмидта. Формула, выражающая турбулентное число Прандтля через характерные времена пульсаций скорости и скаляра, предложена в [37]. Алгебраические модели для турбулентных чисел Прандтля и Шмидта предложены, например, в [38, 39]. Различные дифференциальные модели для турбулентных чисел Прандтля и Шмидта рассмотрены в [40]. Алгебраическая модель [39] и дифференциальная модель [41] предсказывают подобные распределения чисел Прандтля для эксперимента Берроуза и Куркова [42], включающие уменьшение турбулентного числа Прандтля вблизи фронта пламени. Тем не менее, более современная модель [43] не описывает данный эффект. Возможно, ошибкой модели [43] является отождествление характерных частот пульсаций скорости и скалярных параметров. Перспективные модели различных членов в уравнениях, которые используются при расчете турбулентных чисел Прандтля и Шмидта, описаны в работах [44–47].

В работе [48] описан подход к численному моделированию турбулентного горения в струйных течениях, в котором учитываются оба основных канала TCI. Для учета 1-го канала TCI задаются ФПВ концентраций и температуры, а для учета 2-го канала используется модель для переменных турбулентных чисел Прандтля и Шмидта, восходящая к работе [37]. К сожалению, использованный в [48] статистический метод учета 1-го канала TCI не учитывает физические особенности реальных турбулентных пламен; к тому же в [48] не учитываются химические источники в уравнениях для нахождения турбулентных чисел Прандтля и Шмидта.

Современные методы учета TCI требуют значительных компьютерных ресурсов, а точность описания течений с турбулентным горением остается довольно низкой. Разработанные методы описания TCI реализованы в коммерческих пакетах вычислительной аэродинамики (ANSYS CFX [49], FASTRAN [50] и др.), но применение этих методов для решения сложных практических задач требует высокой квалификации вычислителя, глубокого понимания физики течения, умения вторгаться в программу для повышения ее устойчивости и качества.

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

По количеству компонентов и реакций модели химической кинетики можно условно разделить на четыре категории: детальные, скелетные, упрощенные и глобальные. Детальные механизмы горения водорода в воздухе, не учитывающие окисление азота и реакции с участием атомов углерода, обычно содержат около 20 реакций между 9 компонентами (напр., [51]). Детальные механизмы для углеводородов (напр., GriMech 3.0 [52]) включают сотни компонент и тысячи реакций. Скелетные механизмы [53] получаются из детальных путем исключения несущественных для условий задачи веществ и реакций; для углеводородов они могут включать десятки компонент и десятки или сотни реакций. Например, скелетный механизм окисления метана [54] содержит 35 реакций между 15 активными компонентами, а его уточненный аналог [55] – 42 реакций между 17 активными компонентами. Дальнейшее упрощение механизма (за счет сужения области применимости) обеспечивают редуцированные кинетические схемы [56, 57], которые получаются из скелетных с помощью методов квазистационарного состояния и квазиравновесных реакций. Самыми компактными (но и наиболее узкими по применимости) являются глобальные кинетические схемы, которые не выводятся напрямую из других механизмов, а представляют собой небольшой набор брутто-реакций, заменяющих целые этапы реального кинетического процесса. Скорости этих брутто-реакций могут аппроксимироваться абстрактными выражениями, отличающимися от аррениусовых зависимостей (см., напр., [58–60]). Рассматриваются также квазиглобальные модели, в которых одна или несколько глобальных реакций (напр., реакция разложения углеводорода на CO и H₂O) дополнены элементарными реакциями между более простыми частицами (пример – [61]).

Модели химической кинетики играют важную роль в получении точного результата моделирования турбулентного горения. Одновременно с развитием технологий компьютерного моделирования совершались отдельные попытки исследования влияния различных моделей химической кинетики на результат моделирования турбулентного горения. В работе [62] сравнивалось влияние моделей химической кинетики с 8 реакциями и с 25 реакциями на результаты расчетов кругового сверхзвукового потока всасываемого воздуха и двумерной щели в стенке аэродинамической трубы. В работах [63, 64] изучалось влияние неопределенности коэффициентов скорости на вычисленные времена задержки воспламенения, скорости горения и интегральные характеристики течения. Позже в работе [65] были рассмотрены и сопоставлены три кинетики для описания процесса воспламенения неперемешанной системы водород/воздух в сверхзвуковом слое смешения с температурами более 1000 K. В работе [66] на базе осредненных по Рейнольдсу уравнений Навье–Стокса (RANS) исследовано влияние модели химической кинетики на приподнятые сверхзвуковые пламена. Также в работе [67] было исследовано влияние модели химической кинетики на распределение компонентов в пристеночном высокоскоростном горении водорода в воздухе. В последнее десятилетие большое внимание привлекло сочетание метода крупных вихрей (LES) с химическими кинетическими моделями и моделью турбулентного горения. Влияние кинетики на моделирование течения с высокоскоростным горением водорода в канале было рассмотрено в [68–70].

В разд. 2 настоящей статьи представлена базовая система уравнений и замыкающие ее модели турбулентности и химической кинетики. В разд. 3 рассматриваются различные модели 1-го канала TCI класса PaSR. Разд. 4 посвящен модели для переменных турбулентных чисел Прандтля и Шмидта PrOm и ее сопряжению с моделью EPaSR.

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

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

$\frac{\partial \overline{\rho }}{\partial t}+\frac{\partial }{\partial x_i}\left(\overline{\mathrm{\rho }}{\tilde{u}}_i\right)=0$ (2.1)

$\frac{\partial }{\partial t}\left(\overline{\mathrm{\rho }}{\tilde{u}}_k\right)+\frac{\partial }{\partial x_i}\left(\overline{\mathrm{\rho }}{\tilde{u}}_i{\tilde{u}}_k+\overline{p}{\mathrm{\delta }}_{ik}+{\mathrm{\tau }}_{ik}\right)=0$ (2.2)

$\frac{\partial }{\partial t}\left(\overline{\mathrm{\rho }}\tilde{E}\right)+\frac{\partial }{\partial x_i}\left(\overline{\mathrm{\rho }}\tilde{H}{\tilde{u}}_i+{\mathrm{\tau }}_{ik}{\tilde{u}}_k+q_i+\sum _{j=1}^{{\text{N}}_{sp}}{J}_{ij}\mathrm{\Delta }{h}_{\mathrm{0,}j}\right)=0$ (2.3)

$\frac{\partial }{\partial t}\left(\overline{\mathrm{\rho }}{\tilde{Y}}_j\right)+\frac{\partial }{\partial x_i}\left(\overline{\mathrm{\rho }}{\tilde{Y}}_j{\tilde{u}}_i+J_{ij}\right)=\overline{{\dot{s}}_j}$ (2.4)

Здесь t – время, $x_i(i=\mathrm{1,2,3})$ – декартовы координаты, ρ – плотность, $u_i(i=\mathrm{1,2,3})$ – компоненты вектора скорости $\overrightarrow{V}$, $\mathrm{p}=\mathrm{\rho }{R}_{0}T{\sum }_j{Y}_j/W_j$ – давление (R₀ – универсальная газовая постоянная, Y_j и W_j – массовая доля и молекулярный вес j-й компоненты газовой смеси), $E,H=E+p/\mathrm{\rho }$ – полная энергия и полная энтальпия потока на единицу массы газа, ${\dot{s}}_j$ – химический источник массы j-й компоненты, связанный со скоростями химических реакций.

В статье принято суммирование по повторяющимся пространственным индексам. Черта сверху обозначает осреднение по Рейнольдсу, тильда – осреднение по Фавру.

Газ рассматривается как смесь N_sp компонент, так что необходимо решать N_sp уравнений переноса (2.4) для массовых долей ${\tilde{Y}}_i$. Фактически, можно было бы решать только $N_{sp}-1$ уравнений, т.к. сумма ${\tilde{Y}}_i$ всегда равна единице; но ради универсальности алгоритма лучше не исключать уравнение для какой-либо массовой доли. Этот подход требует процедуру перенормировки для поддержания единичной суммы ${\tilde{Y}}_i$ на каждом шаге по времени в течение расчета.

Уравнения (2.2) и (2.3) содержат сумму напряжений Рейнольдса и вязких напряжений, которая представляется как

${\mathrm{\tau }}_{ik}=\overline{\mathrm{\rho }}\stackrel{~}{{u^{\text{'}\text{'}}}_i{u^{\text{'}\text{'}}}_k}-2\tilde{\mathrm{\mu }}\left({\tilde{S}}_{ik}-\frac{1}{3}\frac{\partial {\tilde{u}}_m}{\partial x_m}{\delta }_{ik}\right),$ (2.5)

где ${\tilde{S}}_{ik}=\left(\partial {\tilde{u}}_i/\partial x_k+\partial {\tilde{u}}_k/\partial x_i\right)/2$ – тензор скоростей деформации, m – динамическая молекулярная вязкость, определяемая по формуле Сазерленда. Для воздуха

$\mathrm{\mu }=1.72\cdot {10}^{-5}{\left(\frac{T}{273}\right)}^{3/2}\frac{273+122}{T+122} \frac{\mathrm{кг}}{\mathrm{м}\cdot \mathrm{с}}$

Уравнение (2.3) включает полную энтальпию $\tilde{H}={\tilde{u}}_i{\tilde{u}}_i/2+\tilde{k}+\widetilde{\mathrm{\Delta }{h}_0}+{\tilde{h}}_s$, где $\tilde{k}=\widetilde{u_i^{\text{'}\text{'}}{u}_i^{\text{'}\text{'}}}/2$ – кинетическая энергия турбулентности, $\widetilde{\Delta h_0}=\sum _j{\tilde{Y}}_j\Delta h_{\mathrm{0,}j}$ – удельная энтальпия образования смеси, ${\tilde{h}}_s=\sum _j{\tilde{Y}}_j{h}_{s,j}\left(\tilde{T}\right)$ – явная, или же тепловая энтальпия смеси. Для аппроксимации суммарных энтальпий отдельных компонент, $h_j\left(T\right)=\Delta h_{\mathrm{0,}j}+h_{s,j}\left(T\right)=\Delta h_{\mathrm{0,}j}+\underset{T_0}{\overset{T}{\int }}{c}_{p,j}\left(T^{\text{'}}\right)d{T}^{\text{'}}$, используется база данных [71]. Полная энергия $\tilde{E}$ равна $\tilde{H}-\overline{p}/\overline{\mathrm{\rho }}$. Включение энтальпии образования компонент в $\tilde{H}$ и $\tilde{E}$ позволяет устранить химический источник в правой части уравнения энергии (2.3). Теплоемкость смеси при постоянном давлении считается как $c_p=\sum _j{\tilde{Y}}_{j}d{h}_j/d\tilde{T}$. Также (2.3) содержит сумму турбулентных и молекулярных потоков явной энтальпии:

$q_i=\overline{\mathrm{\rho }}\stackrel{~}{{h^{\text{'}\text{'}}}_s{u^{\text{'}\text{'}}}_i}-\frac{\tilde{\mathrm{\mu }}}{\text{Pr}}\frac{\partial {\tilde{h}}_s}{\partial x_i}$ (2.6)

Уравнения (2.3) и (2.4) содержат также сумму турбулентного и молекулярного потоков массы j-го компонента:

$J_{ji}=\overline{\mathrm{\rho }}\stackrel{~}{{Y^{\text{'}\text{'}}}_j{u^{\text{'}\text{'}}}_i}-\frac{\tilde{\mathrm{\mu }}}{\text{Sc}}\frac{\partial {\tilde{Y}}_j}{\partial x_i}$ (2.7)

Для замыкания системы уравнений (2.1)–(2.4) необходимы три физические модели: модель турбулентности, модель химической кинетики и модель взаимодействия турбулентности и горения (TCI – turbulence combustion interaction).

В рамках данной работы используются дифференциальные модели турбулентности класса $k-\mathrm{\omega }$ (модель SST [72] и модель Baseline $k-\mathrm{\omega }$ [73]). Эти модели основаны на гипотезе Буссинеска, а турбулентные потоки в (2.5)–(2.7) представлены по аналогии с молекулярными потоками:

$\overline{\mathrm{\rho }}\stackrel{~}{{u^{\text{'}\text{'}}}_i{u^{\text{'}\text{'}}}_k}\approx \frac{2}{3}\mathrm{\rho }\tilde{k}{\delta }_{ik}-2{\mathrm{\mu }}_T\left({\tilde{S}}_{ik}-\frac{1}{3}\frac{\partial {\tilde{u}}_m}{\partial x_m}{\mathrm{\delta }}_{ik}\right)$ (2.8)

$\overline{\mathrm{\rho }}\stackrel{~}{{h^{\text{'}\text{'}}}_s{u^{\text{'}\text{'}}}_i}\approx -\frac{{\mathrm{\mu }}_T}{{\text{Pr}}_T}\frac{\partial {\tilde{h}}_s}{\partial x_i}$ (2.9)

$\overline{\mathrm{\rho }}\stackrel{~}{{Y_j}^{\text{'}\text{}\text{'}}{u^{\text{'}\text{'}}}_i}\approx -\frac{{\mathrm{\mu }}_T}{{\text{Sc}}_T}\frac{\partial {\tilde{Y}}_j}{\partial x_i}$ (2.10)

В моделях класса $k-\mathrm{\omega }$ турбулентная вязкость считается следующим образом:

${\mathrm{\mu }}_T=\mathrm{\rho }{F}_{\mathrm{\mu }}\frac{\tilde{k}}{\mathrm{\omega }},$ (2.11)

где ω – характерная частота пульсаций скорости, а $F_{\mathrm{\mu }}$ – пристеночная функция, которая равна единице в свободной турбулентности. В модели SST, $F_{\mathrm{\mu }}=\min \left(1; 0.31\mathrm{\omega }/\left(S{F}_2\right)\right),$ где $S=\sqrt{2S_{ij}{S}_{ij}}$, и $F_2$ – гладкая функция, которая равна единице у стенок и стремится к нулю при удалении от них. В модели Baseline $k-\mathrm{\omega }$, $F_{\mathrm{\mu }}=1$.

Параметры $\tilde{k}$ и ω определяются из решения дополнительных дифференциальных уравнений:

$\begin{array}{l}\frac{\partial }{\partial t}\left(\overline{\mathrm{\rho }}\tilde{k}\right)+\frac{\partial }{\partial x_i}\left(\overline{\mathrm{\rho }}\tilde{k}{\tilde{u}}_i-\left(\tilde{\mathrm{\mu }}+\frac{{\mathrm{\mu }}_T}{{\mathrm{Pr}}_T^k}\right)\frac{\partial \tilde{k}}{\partial x_i}\right)=\overline{\mathrm{\rho }}\left(\tilde{P}-\tilde{\mathrm{\epsilon }}\right)\\ \frac{\partial }{\partial t}\left(\overline{\mathrm{\rho }}\mathrm{\omega }\right)+\frac{\partial }{\partial x_i}\left(\overline{\mathrm{\rho }}\mathrm{\omega }{\tilde{u}}_i-\left(\tilde{\mathrm{\mu }}+\frac{{\mathrm{\mu }}_T}{{\mathrm{Pr}}_T^{\mathrm{\omega }}}\right)\frac{\partial \mathrm{\omega }}{\partial x_i}\right)=\\ = \overline{\mathrm{\rho }}\left(C_{\mathrm{\omega }1}\frac{\tilde{P}}{{\mathrm{\mu }}_T}-C_{\mathrm{\omega }2}\frac{\mathrm{\omega }}{\tilde{k}}\tilde{\epsilon }\right)+C_{\mathrm{\omega }3}\overline{\rho }\left(1-F_1\right)\frac{\partial \tilde{k}}{\partial x_i}\frac{\partial \mathrm{\omega }}{\partial x_i} \end{array}$ (2.12)

$\begin{array}{l}\frac{1}{{\mathrm{Pr}}_T^k}=\frac{1}{{\mathrm{Pr}}_{T\mathrm{,1}}^k}{F}_1+\frac{1}{{\mathrm{Pr}}_{T\mathrm{,2}}^k}(1-F_1),\frac{1}{{\mathrm{Pr}}_T^{\omega }}=\frac{1}{{\mathrm{Pr}}_{T\mathrm{,1}}^{\omega }}{F}_1+\frac{1}{{\mathrm{Pr}}_{T\mathrm{,2}}^{\omega }}(1-F_1)\\ C_{\omega 2}=C_{\omega \mathrm{2,1}}{F}_1+C_{\omega \mathrm{2,2}}(1-F_1),C_{\omega 1}=\frac{C_{\omega 2}}{C_{\mu }}-\frac{{0.41}^2}{{\mathrm{Pr}}_T^{\omega }\sqrt{C_{\mu }}},\end{array}$ (2.13)

где $\tilde{P}=-\stackrel{~}{{u^{\text{'}\text{'}}}_i{u^{\text{'}\text{'}}}_k}\cdot \partial {\tilde{u}}_i/\partial x_k$ – производство кинетической энергии турбулентности $\tilde{k},$ $\tilde{\epsilon }=0.09\tilde{k}\omega$ – скорость диссипации $\tilde{k}$ и $F_1$ – гладкая функция, которая равна единице у стенок и стремится к нулю при удалении от них, ${\mathrm{Pr}}_{T\mathrm{,1}}^k,{\mathrm{Pr}}_{T\mathrm{,1}}^{\omega },C_{\omega \mathrm{2,1}}$ – значения эмпирических коэффициентов для пристеночной турбулентности, ${\mathrm{Pr}}_{T\mathrm{,2}}^k,{\mathrm{Pr}}_{T\mathrm{,2}}^{\omega },C_{\omega \mathrm{2,2}}$ – для свободной турбулентности. Эмпирические коэффициенты ${\mathrm{Pr}}_{T\mathrm{,1}}^k,{\mathrm{Pr}}_{T\mathrm{,1}}^{\omega },C_{\omega \mathrm{2,1}},$ ${\mathrm{Pr}}_{T\mathrm{,2}}^k,{\mathrm{Pr}}_{T\mathrm{,2}}^{\omega },C_{\omega \mathrm{2,2}}$ и $C_{\omega 3}$ указаны в работе [72] (для модели SST) и в работе [73] (для модели Baseline $k-\mathrm{\omega }$). Они отличаются только значениями коэффициента ${\mathrm{Pr}}_{T\mathrm{,1}}^k$, равного 1.176 в модели SST и 2 в Baseline $k-\mathrm{\omega }$.

Модель химической кинетики определяет список веществ и формулы для скоростей реакций, ${\dot{\omega }}_j$, входящих в массовые источники компонентов $\overline{{\dot{s}}_i}$ в уравнении (2.4).

В настоящей работе используются модели обоих основных каналов TCI. Модели 1-го канала TCI (влияние турбулентности на средние скорости химических реакций, ${\overline{\dot{\omega }}}_j$) будут обсуждаться в разд. 3. Если модели 1-го канала TCI не используются, ${\overline{\dot{\omega }}}_j$ можно рассчитать путем подстановки средних параметров в стандартные формулы химических реакций (так называемое квазиламинарное приближение). В рамках гипотезы Буссинеска второй канал TCI можно описать введением переменных турбулентных чисел Прандтля и Шмидта (см. (2.9), (2.10)). Модель для переменных ${\text{Pr}}_T$ и ${\text{Sc}}_T$ будет рассмотрена в разделе 4. Если 2-й канал TCI не учитывается, то расчеты можно выполнить с постоянными числами ${\text{Pr}}_T=\mathrm{0.9,} {\text{Sc}}_T=1.0$.

3. Модели 1-го канала TCI. На основе успешного опыта [24, 25] в данной работе для описания первого канала TCI используется класс моделей реакторов с частичным перемешиванием (PaSR).

Концепция PaSR рассматривает турбулентное течение с горением как двухфазный поток, состоящий из так называемых тонких структур, в которых происходит горение, и окружающей среды, в которой реакции отсутствуют. Экспериментальные данные подтверждают существование тонких структур, см. [20]. Тонкие структуры рассматриваются как реакторы, осуществляющие обмен массой и теплом с окружающей средой за счет диффузии. Параметры в тонких структурах будем обозначать верхним индексом “*”, а параметры окружающей среды индексом “0”. Практичнее решать уравнения для тонких структур и для средних параметров, обозначаемых чертой (осреднение по Рейнольдсу) и тильдой (осреднение по Фавру). Зная объемную долю тонких структур ${\gamma }^{*}$, можно выписать следующие соотношения между этими тремя типами условно осредненных параметров:

$\begin{array}{l}\overline{\mathrm{\rho }}={\gamma }^{*}{\mathrm{\rho }}^{*}+\left(1-{\gamma }^{*}\right){\mathrm{\rho }}^0, \overline{\rho }\widetilde{Y_j}={\gamma }^{*}{\mathrm{\rho }}^{*}{Y}_j^{*}+\left(1-{\gamma }^{*}\right){\mathrm{\rho }}^0{Y}_j^0\\ \overline{\mathrm{\rho }}\sum _{j=1}^{N_{sp}}{\tilde{Y}}_j{h}_j\left(\tilde{T}\right)={\gamma }^{*}{\mathrm{\rho }}^{*}\sum _{j=1}^{N_{sp}}{Y}_j^{*}{h}_j\left(T^{*}\right)+\left(1-{\gamma }^{*}\right){\mathrm{\rho }}^0\sum _{j=1}^{N_{sp}}{Y}_j^0{h}_j\left(T^0\right)\end{array}$ (3.1)

Наиболее известными моделями этого класса являются Eddy Dissipation Concept (EDC) [20, 74] и PaSR [21]. Эти модели основаны на предположении, что время реакции в тонких структурах существенно меньше характерного времени изменения среднего потока. В этом случае состояние тонких структур является квазистационарным и описывается следующей системой алгебраических уравнений:

$\begin{array}{l}{\mathrm{\rho }}^{*}\frac{Y_j^{*}-Y_j^{\text{0}}}{{\tau }^{*}}={\dot{s}}_j\left({\mathrm{\rho }}^{*},T^{*},Y_1^{*},\dots ,Y_{N_{sp}}^{*}\right)\\ \sum _{j=1}^{N_{sp}}{Y}_j^{*}{h}_j\left(T^{*}\right)=\sum _{j=1}^{N_{sp}}{Y}_j^0{h}_j\left(T^0\right)\end{array}$ (3.2)

Значение ${\mathrm{\rho }}^{*}/{\mathrm{\tau }}^{*}$ является средним потоком массы через единичный объем тонких структур, и ${\mathrm{\tau }}^{*}$ может быть интерпретировано как характерное время смешения, определяющее скорость обмена тонких структур массой и теплом с окружающим пространством. Скорость газа, как в тонких структурах, так и в окружающей среде предполагается равной локальной средней скорости потока. Следовательно, естественно принять, что

$p^{*}=p^0=\overline{p}$ (3.3)

Соотношение (3.3) используется в реализации подхода PaSR в программе zFlare, разработанной в ЦАГИ [75]. Тем не менее, в программе ANSYS Fluent© [49] используется другое предположение вместо (3.3):

${\mathrm{\rho }}^{*}={\mathrm{\rho }}^0=\overline{\mathrm{\rho }}$ (3.4)

Уравнения (3.1)–(3.2) образуют замкнутую систему алгебраических уравнений, позволяющую через заданные параметры осредненного потока $\left(\overline{\mathrm{\rho }},\tilde{T},{\tilde{Y}}_1\mathrm{,...}{\tilde{Y}}_{N_{sp}}\right)$ найти параметры газа в тонких структурах и в окружающей среде. После нахождения параметров газа в тонких структурах средние скорости химических реакций определяются по следующей формуле:

$\overline{{\dot{s}}_j}\approx {\gamma }^{*}{\dot{s}}_j\left({\mathrm{\rho }}^{*},T^{*},Y_1^{*},\dots ,Y_{N_{sp}}^{*}\right)={\gamma }^{*}{\mathrm{\rho }}^{*}\frac{Y_j^{*}-Y_j^{\text{0}}}{{\tau }^{*}}$ (3.5)

В модели PaSR [21] объемная доля тонких структур ${\gamma }^{*}$ и характерное время смешения газа в тонких структурах выражаются формулами

${\mathrm{\gamma }}^{*}=\frac{{\mathrm{\tau }}_{\mathrm{chem}}}{{\mathrm{\tau }}_{\mathrm{chem}}+{\mathrm{\tau }}^{*}}, {\mathrm{\tau }}^{*}={\mathrm{C}}_{\mathrm{\tau }}\sqrt{{\mathrm{\tau }}_T{\mathrm{\tau }}_K}$, (3.6)

где ${\mathrm{\tau }}_{\mathrm{chem}}$ – характерное время реакции в тонких структурах, ${\mathrm{\tau }}_T=\tilde{k}/\tilde{\mathrm{\epsilon }}~1/\mathrm{\omega }$ – интегральный временной масштаб турбулентных пульсаций, а ${\mathrm{\tau }}_K=\sqrt{\tilde{\nu }/\tilde{\mathrm{\epsilon }}}$ – микромасштаб Колмогорова. Отметим, что формулу (3.6) можно переписать в виде ${\mathrm{\gamma }}^{*}={\left[1+C_{\mathrm{\tau }}\sqrt{\mathrm{Da}/\mathrm{Ka}}\right]}^{-1}$, где $\mathrm{Da}={\mathrm{\tau }}_T/{\mathrm{\tau }}_{\mathrm{chem}}$ и $\mathrm{Ka}={\mathrm{\tau }}_{\mathrm{chem}}/{\mathrm{\tau }}_K=\sqrt{{\mathrm{Re}}_T}/\mathrm{Da}$ – числа Дамкелера и Карловица, сочетание которых определяет режим турбулентного горения [4]. Коэффициент $C_{\mathrm{\tau }}$в работе [21] принят равным единице.

Решение системы алгебраических уравнений (3.1)–(3.2) в zFlare реализовано с использованием установления по псевдовремени. Следующие уравнения интегрируются по псевдовремени $t^{\text{'}}$ до достижения стационарного состояния ($d/d{t}^{\text{'}}\to 0$):

$\begin{array}{l}{\mathrm{\rho }}^{*}\frac{d{Y}_j^{*}}{d{t}^{\text{'}}}+{\mathrm{\rho }}^{*}\frac{Y_j^{*}-Y_j^0}{{\tau }^{*}}={\dot{s}}_j\left({\mathrm{\rho }}^{*},T^{*},Y_1^{*},\dots ,Y_{N_{sp}}^{*}\right)\\ \sum _{j=1}^{N_{sp}}{Y}_j^{*}{h}_j\left(T^{*}\right)=\sum _{j=1}^{N_{sp}}{Y}_j^0{h}_j\left(T^0\right)\end{array}$ (3.7)

При интегрировании по псевдовремени используется формула (3.3).

В ANSYS Fluent© используется модификация модели PaSR – так называемая PFR (Plug Flow Reactor [76–78]), где (3.2) с допущением (3.4) заменяется на систему обыкновенных дифференциальных уравнений:

$\begin{array}{l}\overline{\mathrm{\rho }}\frac{d{Y}_j^{*}}{d{t}^{\text{'}}}={\dot{s}}_j\left(\overline{\text{ρ}},T^{*},Y_1^{*},\dots ,Y_{N_{sp}}^{*}\right)\\ \sum _{j=1}^{N_{sp}}{Y}_j^{*}{h}_j\left(T^{*}\right)=\sum _{j=1}^{N_{sp}}{Y}_j^0{h}_j\left(T^0\right)\end{array}$ (3.8)

В начальный момент псевдовремени, $t^{\text{'}}=0$, параметры в тонких структурах принимаются равными параметрам окружающего газа: $Y_j^{*}\left(0\right)=Y_j^0$, $T^{*}\left(0\right)=T^0$. Интегрирование (3.8) по псевдовремени ведется в интервале $t^{\text{'}}\in \left[0; {\tau }^{*}\right]$. В итоге, вместо (3.5) используется следующая формула:

$\overline{{\dot{s}}_j}\approx {\gamma }^{*}\underset{0}{\overset{{\tau }^{*}}{\int }}{\dot{s}}_j\left(\overline{\mathrm{\rho }},T^{*}\left(t^{\text{'}}\right),Y_1^{*}\left(t^{\text{'}}\right),\dots ,Y_{{\text{N}}_{sp}}^{*}\left(t^{\text{'}}\right)\right)d{t}^{\text{'}}={\gamma }^{*}\overline{\mathrm{\rho }}\frac{Y_j^{*}\left({\tau }^{*}\right)-Y_j^{\text{0}}}{{\tau }^{*}}$ (3.9)

В работах [22, 79] была предложена расширенная модель PaSR (Extended PaSR – EPaSR). Для учета изменения параметров тонких структур в пространстве и во времени было предложено решать дополнительные уравнения в частных производных для параметров тонких структур:

$\begin{array}{l}\frac{\partial \left({\gamma }^{*}{\mathrm{\rho }}^{*}{Y}_j^{*}\right)}{\partial t}+\frac{\partial }{\partial x_i}\left[{\gamma }^{*}\left({\mathrm{\rho }}^{*}{Y}_j^{*}{\tilde{u}}_i+J_{ji}^{*}\right)\right]=\\ =-{\gamma }^{*}{\mathrm{\rho }}^{*}\left(C_g\frac{\tilde{\epsilon }}{\tilde{k}}\left(Y_j^{*}-Y_j^0\right)\right)+{\gamma }^{*}{\mathrm{\rho }}^{*}\overline{{\dot{s}}_j}\left({\mathrm{\rho }}^{*},T^{*},Y_1^{*},\dots ,Y_{N_{sp}}^{*}\right)+\\ +\frac{\dot{m}+\left|\dot{m}\right|}{2}{Y}_j^0+\frac{\dot{m}-\left|\dot{m}\right|}{2}{Y}_j^{*};j=\mathrm{1,...,}{N}_{sp}-1\end{array}$ (3.10)

$\begin{array}{l}\frac{\partial \left({\gamma }^{*}{\mathrm{\rho }}^{*}{E}^{*}\right)}{\partial t}+\frac{\partial }{\partial x_i}\left[{\gamma }^{*}\left({\mathrm{\rho }}^{*}{E}^{*}{\tilde{u}}_i+\overline{p}{\tilde{u}}_i+{\tau }_{ki}{\tilde{u}}_k+q_i^{*}+\sum _{i=1}^{N_{sp}}{J}_{ji}^{*}\Delta h_{\mathrm{0,}j})\right)\right]=\\ =-{\gamma }^{*}{\mathrm{\rho }}^{*}{C}_g\frac{\tilde{\epsilon }}{\tilde{k}}\sum _{j=1}^{N_{sp}}(Y_j^{*}{h}_j\left(T^{*}\right)-Y_j^0{h}_j\left(T^0\right))+\\ +\frac{\dot{m}+\left|\dot{m}\right|}{2}\sum _{j=1}^{N_{sp}}{Y}_j^0{h}_j\left(T^0\right)+\frac{\dot{m}-\left|\dot{m}\right|}{2}\sum _{j=1}^{N_{sp}}{Y}_j^{*}{h}_j\left(T^{*}\right)\end{array}$ (3.11)

$\frac{\partial \left({\gamma }^{*}{\mathrm{\rho }}^{*}\right)}{\partial t}+\frac{\partial }{\partial x_i}\left[{\gamma }^{*}{\mathrm{\rho }}^{*}{\tilde{u}}_i\right]=\dot{m}; \dot{m}=-{\mathrm{\rho }}^{*}\frac{{\gamma }^{*}-{\gamma }_{eq}^{*}}{{\tau }^{*}}$ (3.12)

Константа $C_g$ в релаксационных терминах (описывающая диффузионный обмен между тонкими структурами и окружающей средой) была выбрана равной 10.5, исходя из теоретических соображений [24]. Это было подтверждено параметрическими расчетами [80].

Слагаемые с множителем $\dot{m}$ ответственны за изменение параметров вследствие обмена массы от тонких структур к окружающему газу и обратно. Объемная доля тонких структур, ${\gamma }^{*}$, определяется из уравнения (3.12), где ${\mathrm{\gamma }}_{\mathrm{eq}}^{*}={\mathrm{\tau }}_{\text{chem}}/\left({\mathrm{\tau }}_{\text{chem}}+{\mathrm{\tau }}^{*}\right)$. На практике вместо решения уравнения (3.12) уравнение (3.10) решается также для $j=N_{sp}$. (Сумма уравнений (3.10) для $j=\mathrm{1,2,...,}{N}_{sp}$ эквивалентна уравнению (3.12).) Вместе с системой (3.10)–(3.12) используется соотношение (3.3). Также предполагается, что турбулентная вязкость и турбулентные числа Прандтля и Шмидта для тонких структур такие же, как и для среднего течения, и используются следующие формулы для диффузионных потоков:

$q_i^{*}=-\left(\frac{\mathrm{\mu }\left(T^{*}\right)}{\text{Pr}}+\frac{{\mathrm{\mu }}_T}{{\text{Pr}}_T}\right)\frac{\partial }{\partial x_i}\left(\sum _{j=1}^{N_{sp}}{Y}_j^{*}{h}_{s,j}\left(T^{*}\right)\right)$ (3.13)

$J_{ji}^{*}=-\left(\frac{\mathrm{\mu }\left(T^{*}\right)}{\text{Sc}}+\frac{{\mathrm{\mu }}_T}{{\text{Sc}}_T}\right)\frac{\partial Y_j^{*}}{\partial x_i}$ (3.14)

4. Модели 2-го канала TCI. Модель PrOm [81] описывает турбулентный перенос на основе переменных турбулентных чисел Прандтля и Шмидта для потоков тепла и массы соответственно. Этот подход рассматривает турбулентное число Прандтля для произвольного скаляра f как функцию частоты пульсаций скорости ω и частоты пульсаций скаляра ω_f. Первая частота находится из уравнений для модели турбулентности класса k - ω ( (Baseline k - ω или SST). Для нахождения частоты пульсаций скаляра ω_f требуется решение дополнительных дифференциальных уравнений для среднего квадрата пульсаций скаляра $\stackrel{~}{{f^{\text{'}\text{'}}}^2}$, а также для самой частоты пульсаций скаляра ω_f. Более предпочтительно решать отдельные дифференциальные уравнения для частот, поскольку асимптотическое поведение их решения вблизи стенки близко к реальному [2], и их низкорейнольдсовые версии не требуют искусственных эмпирических пристеночных коррекций.

Формула для турбулентного числа Прандтля была впервые предложена в работе [37]. Существуют различные способы получения этой формулы. Один из способов представлен в Приложении. Он основан на рассмотрении дифференциального уравнения для турбулентного потока скаляра $\stackrel{~}{f^{\text{'}\text{'}}{u^{\text{'}\text{'}}}_i}$. Обзор различных замыканий для этого уравнения представлен в Приложении. В результате применения этих замыканий и приравнивания нулю правой части уравнения для $\stackrel{~}{f^{\text{'}\text{'}}{u^{\text{'}\text{'}}}_i}$ можно получить следующую простую формулу для турбулентного числа Прандтля в приближении Буссинеска $\stackrel{~}{f^{\text{'}\text{'}}{u^{\text{'}\text{'}}}_i}\approx -\frac{{\mathrm{\mu }}_T}{\overline{\mathrm{\rho }}{\mathrm{Pr}}_T^f}\frac{\partial \tilde{f}}{\partial x_i}$:

$P{r}_T^f=C_{\Phi }\sqrt{\frac{{\omega }_f}{\omega }}$ (4.1)

Для простоты определим частоту пульсаций скаляра по формуле, аналогичной формуле для частоты пульсаций скорости $\mathrm{\omega }=\tilde{\mathrm{\epsilon }}/\left(C_{\mathrm{\mu }}\tilde{k}\right),C_{\mathrm{\mu }}=0.09$:

${\mathrm{\omega }}_f=\frac{{\mathrm{\epsilon }}_f}{0.09{\mathrm{\sigma }}_f},$ (4.2)

где ${\mathrm{\sigma }}_f\equiv \stackrel{~}{{f^{\text{'}\text{'}}}^2}$ – квадрат пульсаций скаляра и ${\epsilon }_f=2\frac{\tilde{\mathrm{\mu }}}{\overline{\mathrm{\rho }}}\overline{\frac{\partial f^{\text{'}\text{'}}}{\partial x_j}\frac{\partial f^{\text{'}\text{'}}}{\partial x_j}}$ – его диссипация.

Согласно измерениям [82] соотношение частот пульсаций скаляра и скорости варьируется в диапазоне от 0.67 до 2.38. По данным Сринивасана и др., в экспериментах для сдвиговых течений [83] число Прандтля варьируется от 0.4 до 1.1. Предположим, что наименьшее турбулентное число Прандтля 0.4 достигается при ${\mathrm{\omega }}_f=0.67\mathrm{\omega }$, тогда, согласно формуле (4.1), $C_{\mathrm{\Phi }}=0.51$. Это значение будет в дальнейшем использоваться в данной статье.

Стоит отметить, что существует сильный разброс в оценке этой константы, например, в работе [37] предлагается значение $C_{\mathrm{\Phi }}$, равное 0.82. В диссертации [48], значение $C_{\mathrm{\Phi }}=0.52$ было получено путем применения соотношения между квадратом пульсации нормальной к линии тока компоненты скорости и кинетической энергией турбулентности: $\stackrel{~}{{{v^{\text{'}\text{'}}}^2}_n}/\tilde{k}\approx 0.5$.

Точные уравнения для квадрата пульсации скаляра ${\mathrm{\sigma }}_f$ и ${\mathrm{\epsilon }}_f$ его диссипации приведены ниже:

$\begin{array}{l}\frac{\partial \overline{\mathrm{\rho }}{\mathrm{\sigma }}_f}{\partial t}+\frac{\partial }{\partial x_i}\left[\overline{\mathrm{\rho }}{\mathrm{\sigma }}_f{\tilde{u}}_i+\overline{\mathrm{\rho }}\stackrel{~}{{f^{\text{'}\text{'}}}^2{u^{\text{'}\text{'}}}_i}-\tilde{\mathrm{\mu }}\frac{\partial {\mathrm{\sigma }}_f}{\partial x_i}\right]=\\ =-2\overline{\mathrm{\rho }}\stackrel{~}{f^{\text{'}\text{'}}{u^{\text{'}\text{'}}}_k}\frac{\partial \tilde{f}}{\partial x_k}-2\tilde{\mathrm{\mu }}\overline{\frac{\partial f^{\text{'}\text{'}}}{\partial x_k}\frac{\partial f^{\text{'}\text{'}}}{\partial x_k}}+2\overline{f^{\text{'}\text{'}}{\dot{s}}_f}\end{array}$ (4.3)

$\begin{array}{l}\frac{\partial \overline{\mathrm{\rho }}{\epsilon }_f}{\partial t}+\frac{\partial }{\partial x_i}\left[\overline{\mathrm{\rho }}{\epsilon }_f{\tilde{u}}_i+\overline{\mathrm{\rho }}\stackrel{~}{{{\epsilon }^{\text{'}\text{'}}}_f{u^{\text{'}\text{'}}}_i}-\tilde{\mathrm{\mu }}\frac{\partial {\epsilon }_f}{\partial x_i}\right]=\\ =-\underset{\text{EI}}{\underbrace{4\tilde{\mathrm{\mu }}\overline{\frac{\partial f^{\text{'}\text{'}}}{\partial x_m}\frac{\partial f^{\text{'}\text{'}}}{\partial x_k}}\frac{\partial {\tilde{u}}_k}{\partial x_m}}}-\underset{\text{EII}}{\underbrace{4\tilde{\mathrm{\mu }}\overline{\frac{\partial f^{\text{'}\text{'}}}{\partial x_m}\frac{\partial {u^{\text{'}\text{'}}}_k}{\partial x_m}}\frac{\partial \tilde{f}}{\partial x_k}}}-\end{array}$ (4.4)

$\underset{\text{EIII}}{\underbrace{-4\tilde{\mathrm{\mu }}\cdot \frac{\tilde{\mathrm{\mu }}}{\overline{\mathrm{\rho }}}\overline{\frac{{\partial }^2{f}^{\text{'}\text{'}}}{\partial x_m\partial x_k}\frac{{\partial }^2{f}^{\text{'}\text{'}}}{\partial x_m\partial x_k}}-4\tilde{\mathrm{\mu }}\overline{\frac{\partial f^{\text{'}\text{'}}}{\partial x_m}\frac{\partial f^{\text{'}\text{'}}}{\partial x_k}\frac{\partial {u^{\text{'}\text{'}}}_k}{\partial x_m}}}}+\underset{\text{EIV}}{\underbrace{4\frac{\tilde{\mathrm{\mu }}}{\overline{\mathrm{\rho }}}\overline{\frac{\partial f^{\text{'}\text{'}}}{\partial x_m}\frac{\partial {{\dot{s}}^{\text{'}\text{'}}}_f}{\partial x_m}}}}$

В правой части уравнения для скалярной диссипации (4.4) имеются четыре незамкнутых члена. Слагаемое $\text{EI}$ ответственно за производство ${\epsilon }_f$ градиентами среднего поля скорости. Однако, это одно из самых противоречивых слагаемых, так как оно может также вести себя как сток, особенно для положительной дивергенции скорости и для расширяющихся топологий тензора градиента скорости $\partial {\tilde{u}}_i/\partial x_j$ (как было показано в [45]). Замыкание, предложенное в текущей статье, разрешает отрицательные значения этого слагаемого при помощи добавления к стандартной аппроксимации (основанной на производстве кинетической энергии турбулентности $\tilde{P}$) дополнительного слагаемого, отрицательного и квадратичного по дивергенции скорости, как показано ниже:

$\text{EI}=\overline{\mathrm{\rho }}\frac{{\epsilon }_f}{\tilde{k}}\left\{C_{P2}\tilde{P}-C_{P3}{\mathrm{\mu }}_T{\left(\frac{\partial {\tilde{u}}_m}{\partial x_m}\right)}^2\right\}$

Здесь $C_{P2}=1.45$ [40]. Значение коэффициента $C_{P3}$ было получено авторами в параметрических расчетах эксперимента [42]. Второе слагаемое препятствует чрезмерному росту ${\text{Pr}}_T^f$ в областях с повышенным значениями дивергенции скорости, например, области с косым скачком уплотнения и отраженной волной. Слагаемое, обозначенное как $\text{EII}$, представляет собой производство ${\epsilon }_f$ вследствие среднего градиента скаляра и является сверткой среднего градиента скаляра с вектором диссипации турбулентного потока скаляра: ${\text{FIII}}_{\text{k}}=4\mathrm{\mu }\overline{\partial f^{\text{'}\text{'}}/\partial x_m\partial {u^{\text{'}\text{'}}}_k/\partial x_m}$, рассмотренным в Приложении. Пренебрегая вторым слагаемым в выражении для ${\text{FIII}}_{\text{k}}$ (П.2), и применяя формулу Буссинеска для потока скаляра, можно предложить следующее замыкание для слагаемого $\text{EII}$:

$\text{EII}=C_{P1}\frac{{\epsilon }_f}{{\mathrm{\sigma }}_f}2\frac{{\mathrm{\mu }}_T}{{\text{Pr}}_T^f}\frac{\partial \tilde{f}}{\partial x_m}\frac{\partial \tilde{f}}{\partial x_m}=C_{P1}\overline{\rho }\frac{{\epsilon }_f}{{\mathrm{\sigma }}_f}{P}_f, C_{P1}=\frac{0.46}{{{\text{(Pr}}_T^f)}^2}$

Два отдельных слагаемых в $\text{EIII}$ представляют собой деструкцию ${\epsilon }_f$ из-за кривизны изоповерхностей поля скаляра и за счет механизма турбулентного растяжения, соответственно. Замыкание этих членов для режима флеймлетов рассматривалось в [47]. Авторы считают, что слагаемое, ответственное за турбулентное растяжение, должно быть линейным по скалярной диссипации, тогда как член, связанный с кривизной изоповерхностей скаляра, должен быть квадратичным. Несмотря на наличие зависимости от турбулентного числа Рейнольдса, эта модель упускает один важный факт, который был получен для пассивного скаляра в [44]. Согласно экспериментальным данным для сдвиговых потоков [83] анизотропия поля средней скорости приводит к резкому уменьшению турбулентного числа Прандтля. С другой стороны, Ши, Ламли и Яника [44] вывели зависимость квадратичной деструкции от второго инварианта $II$ тензора анизотропии $a_{ij}$, где

$a_{ij}=\frac{\stackrel{~}{{u^{\text{'}\text{'}}}_i{u^{\text{'}\text{'}}}_j}}{k}-\frac{2}{3}{\mathrm{\delta }}_{ij}, II=a_{ij}{a}_{ij}$

Принимая во внимание это исследование, предложенное замыкание для члена $\text{EIII}$ имеет вид:

$\text{EIII}=-\overline{\mathrm{\rho }}{\epsilon }_f\left\{C_{D1}\frac{{\epsilon }_f}{{\mathrm{\sigma }}_f}+C_{D2}\overline{\mathrm{\rho }}\frac{\tilde{\epsilon }}{\tilde{k}}\right\}; C_{D1}=1+\frac{3}{4}{b}_{ij}{b}_{ij}, b_{ij}=\frac{a_{ij}}{2}, C_{D2}=0.9$

Тензор анизотропии рассчитывается по формуле Буссинеска (2.8). Следовательно, первый коэффициент в замыкании выражается как

$C_{D1}=1+\frac{3}{4}{S}_{ij}^{dil}{S}_{ij}^{dil}/{\mathrm{\omega }}^2$,

где $S_{ij}^{dil}={\tilde{S}}_{ij}-\frac{1}{3}\frac{\partial {\tilde{u}}_n}{\partial x_n}{\mathrm{\delta }}_{ij}$.

Слагаемое $\text{EIV}$ рассматривалось в [46] для активного скаляра в случае предварительно перемешанного пламени и выражается через скорость ламинарного фронта пламени и включается в деструкцию на основе кривизны изоповерхностей скаляра, квадратичную по скалярной диссипации.

Другой подход, который используется в текущей статье, заключается в том, чтобы выразить $\text{EIV}$ как функцию от химического источника $Q_f^{\mathrm{chem}}\equiv 2\overline{f^{\text{'}\text{'}}{\dot{s}}_f}$. Итоговое замыкание записывается следующим образом:

$\text{EIV}=C_{W1}\frac{{\epsilon }_f}{{\mathrm{\sigma }}_f}{Q}_f^{\mathrm{chem}}; C_{W1}=0.4$

Замыкание $Q_f^{\mathrm{chem}}$ на основе модели EPaSR рассмотрено ниже.

Уравнения для частоты пульсаций скаляра ${\mathrm{\omega }}_f$ могут быть легко выведены из (4.3) и (4.4), применяя соотношение (4.2).

Ради краткости, мы, наконец, рассмотрим уравнения, описывающие турбулентный перенос явной энтальпии $h_s$ и массовых долей $Y_i$. Согласно (2.9) и (2.10), они требуют определения турбулентного числа Прандтля ${\text{Pr}}_T$ и турбулентного числа Шмидта ${\text{Sc}}_T$. Эти числа могут быть посчитаны по формулам, аналогичным (4.1):

${\text{Pr}}_T=C_{\mathrm{\Phi }}\sqrt{\frac{{\mathrm{\omega }}_h}{\mathrm{\omega }}}, {\text{Sc}}_T=C_{\mathrm{\Phi }}\sqrt{\frac{{\mathrm{\omega }}_C}{\mathrm{\omega }}}$ (4.5)

В этих формулах (см. (4.2)):

${\mathrm{\omega }}_h=\frac{{\mathrm{\epsilon }}_h}{0.09{\mathrm{\sigma }}_h}, {\mathrm{\omega }}_C=\frac{{\mathrm{\epsilon }}_C}{0.09{\mathrm{\sigma }}_C}$ (4.6)

Здесь ${\mathrm{\sigma }}_h=\stackrel{~}{{h^{\text{'}\text{'}}}_s{}^2}$, ${\mathrm{\sigma }}_C=\sum _{j=1}^{N_{sp}}\stackrel{~}{{Y^{\text{'}\text{'}}}_j{}^2}$,  и ${\epsilon }_h$, ${\epsilon }_C$ – их диссипации соответственно. Следовательно, определение турбулентного числа Прандтля ${\text{Pr}}_T$ и турбулентного числа Шмидта ${\text{Sc}}_T$ требует решения четырех дополнительных дифференциальных уравнений для ${\mathrm{\sigma }}_h$, ${\mathrm{\omega }}_h$, ${\mathrm{\sigma }}_C$ и ${\mathrm{\omega }}_C$. Эти уравнения приведены ниже:

$\frac{\partial \overline{\mathrm{\rho }}{\mathrm{\sigma }}_{\mathrm{h}}}{\partial \mathrm{t}}+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{k}}}\left[\overline{\mathrm{\rho }}{\mathrm{\sigma }}_{\mathrm{h}}{\tilde{\mathrm{u}}}_{\mathrm{k}}-{\mathrm{D}}_{\mathrm{T}}^{\mathrm{h}}\frac{\partial {\mathrm{\sigma }}_{\mathrm{h}}}{\partial {\mathrm{x}}_{\mathrm{k}}}\right]=2\frac{{\mathrm{\mu }}_{\mathrm{T}}}{{\mathrm{Pr}}_{\mathrm{T}}^{\mathrm{h}}}{\left(\frac{\partial \widetilde{{\mathrm{h}}_{\mathrm{s}}}}{\partial {\mathrm{x}}_{\mathrm{k}}}\right)}^2-\overline{\mathrm{\rho }}{\mathrm{\epsilon }}_{\mathrm{h}}+2\overline{{{\mathrm{h}}^{\text{'}\text{'}}}_{\mathrm{s}}{\dot{\mathrm{s}}}_{\mathrm{h}}}$ (4.7)

$\begin{array}{l}\frac{\partial \overline{\mathrm{\rho }}{\mathrm{\omega }}_{\mathrm{h}}}{\partial \mathrm{t}}+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{k}}}\left[\overline{\mathrm{\rho }}{\mathrm{\omega }}_{\mathrm{h}}{\tilde{\mathrm{u}}}_{\mathrm{k}}-{\mathrm{D}}_{\mathrm{T}}^{\mathrm{h}}\frac{\partial {\mathrm{\omega }}_{\mathrm{h}}}{\partial {\mathrm{x}}_{\mathrm{k}}}\right]=\overline{\mathrm{\rho }}\frac{{\mathrm{\omega }}_{\mathrm{h}}}{{\mathrm{\sigma }}_{\mathrm{h}}}\left({\mathrm{C}}_{\mathrm{P}1}^{\mathrm{h}}-1\right)2\frac{{\mathrm{\mu }}_{\mathrm{T}}}{{\mathrm{Pr}}_{\mathrm{T}}^{\mathrm{h}}}{\left(\frac{\partial \widetilde{{\mathrm{h}}_{\mathrm{s}}}}{\partial {\mathrm{x}}_{\mathrm{k}}}\right)}^2+\\ +{\mathrm{C}}_{\mathrm{P}2}{\mathrm{\mu }}_{\mathrm{T}}\frac{{\mathrm{\omega }}_{\mathrm{h}}}{\tilde{\mathrm{k}}}2{\mathrm{S}}_{\mathrm{ij}}^{\mathrm{dil}}{\mathrm{S}}_{\mathrm{ij}}^{\mathrm{dil}}-{\mathrm{C}}_{\mathrm{P}3}{\mathrm{\mu }}_{\mathrm{T}}\frac{{\mathrm{\omega }}_{\mathrm{h}}}{\tilde{\mathrm{k}}}{\left(\frac{\partial {\tilde{\mathrm{u}}}_{\mathrm{m}}}{\partial {\mathrm{x}}_{\mathrm{m}}}\right)}^2+\\ +\left(1-{\mathrm{C}}_{\mathrm{D}1}\right){\mathrm{C}}_{\mathrm{\mu }}\overline{\mathrm{\rho }}{\mathrm{\omega }}_{\mathrm{h}}^2-{\mathrm{C}}_{\mathrm{D}2}\overline{\mathrm{\rho }}\frac{\tilde{\mathrm{\epsilon }}{\mathrm{\omega }}_{\mathrm{h}}}{\tilde{\mathrm{k}}}-2\left(1-\frac{{\mathrm{C}}_{\mathrm{W}1}}{2}\right)\frac{{\mathrm{\omega }}_{\mathrm{h}}}{{\mathrm{\sigma }}_{\mathrm{h}}}\overline{{{\mathrm{h}}^{\text{'}\text{'}}}_{\mathrm{s}}{\dot{\mathrm{s}}}_{\mathrm{h}}}\end{array}$ (4.8)

$\frac{\partial \overline{\mathrm{\rho }}{\mathrm{\sigma }}_{\mathrm{C}}}{\partial \mathrm{t}}+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{k}}}\left[\overline{\mathrm{\rho }}{\mathrm{\sigma }}_{\mathrm{C}}{\tilde{\mathrm{u}}}_{\mathrm{k}}-{\mathrm{D}}_{\mathrm{T}}^{\mathrm{C}}\frac{\partial {\mathrm{\sigma }}_{\mathrm{C}}}{\partial {\mathrm{x}}_{\mathrm{k}}}\right]=2\frac{{\mathrm{\mu }}_{\mathrm{T}}}{{\mathrm{Sc}}_{\mathrm{T}}}\sum _{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{sp}}}\frac{\partial \widetilde{{\mathrm{Y}}_{\mathrm{j}}}}{\partial {\mathrm{x}}_{\mathrm{k}}}\frac{\partial \widetilde{{\mathrm{Y}}_{\mathrm{j}}}}{\partial {\mathrm{x}}_{\mathrm{k}}}-\overline{\mathrm{\rho }}{\mathrm{\epsilon }}_{\mathrm{C}}+2\sum _{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{sp}}}\overline{{{\mathrm{Y}}^{\text{'}\text{'}}}_{\mathrm{j}}{\dot{\mathrm{s}}}_{\mathrm{j}}}$ (4.9)

$\begin{array}{l}\frac{\partial \overline{\mathrm{\rho }}{\mathrm{\omega }}_{\mathrm{C}}}{\partial \mathrm{t}}+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{k}}}\left[\overline{\mathrm{\rho }}{\mathrm{\omega }}_{\mathrm{C}}{\tilde{\mathrm{u}}}_{\mathrm{k}}-{\mathrm{D}}_{\mathrm{T}}^{\mathrm{C}}\frac{\partial {\mathrm{\omega }}_{\mathrm{h}}}{\partial {\mathrm{x}}_{\mathrm{k}}}\right]=\overline{\mathrm{\rho }}\frac{{\mathrm{\omega }}_{\mathrm{C}}}{{\mathrm{\sigma }}_{\mathrm{C}}}\left({\mathrm{C}}_{\mathrm{P}1}^{\mathrm{C}}-1\right)2\frac{{\mathrm{\mu }}_{\mathrm{T}}}{{\mathrm{Sc}}_{\mathrm{T}}}\sum _{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{sp}}}\frac{\partial \widetilde{{\mathrm{Y}}_{\mathrm{j}}}}{\partial {\mathrm{x}}_{\mathrm{k}}}\frac{\partial \widetilde{{\mathrm{Y}}_{\mathrm{j}}}}{\partial {\mathrm{x}}_{\mathrm{k}}}+\\ +{\mathrm{C}}_{\mathrm{P}2}{\mathrm{\mu }}_{\mathrm{T}}\frac{{\mathrm{\omega }}_{\mathrm{C}}}{\tilde{\mathrm{k}}}2{\mathrm{S}}_{\mathrm{ij}}^{\mathrm{dil}}{\mathrm{S}}_{\mathrm{ij}}^{\mathrm{dil}}-{\mathrm{C}}_{\mathrm{P}3}{\mathrm{\mu }}_{\mathrm{T}}\frac{{\mathrm{\omega }}_{\mathrm{C}}}{\tilde{\mathrm{k}}}{\left(\frac{\partial {\tilde{\mathrm{u}}}_{\mathrm{m}}}{\partial {\mathrm{x}}_{\mathrm{m}}}\right)}^2+\\ +\left(1-{\mathrm{C}}_{\mathrm{D}1}\right)\overline{\mathrm{\rho }}{\mathrm{C}}_{\mathrm{\mu }}{\mathrm{\omega }}_{\mathrm{C}}^2-{\mathrm{C}}_{\mathrm{D}2}\overline{\mathrm{\rho }}\frac{{\mathrm{\epsilon \omega }}_{\mathrm{C}}}{\tilde{\mathrm{k}}}-2\left(1-\frac{{\mathrm{C}}_{\mathrm{W}1}}{2}\right)\frac{{\mathrm{\omega }}_{\mathrm{C}}}{{\mathrm{\sigma }}_{\mathrm{C}}}\sum _{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{sp}}}\overline{{{\mathrm{Y}}^{\text{'}\text{'}}}_{\mathrm{j}}{\dot{\mathrm{s}}}_{\mathrm{j}}}\end{array}$ (4.10)

Здесь использованы обозначения $D_T^h=\tilde{\mathrm{\mu }}/\text{Pr}+{\mathrm{\mu }}_T/{\text{Pr}}_{\text{T}}$ и $D_T^C=\tilde{\mathrm{\mu }}/\text{Sc}+{\mathrm{\mu }}_T/{\text{Sc}}_{\text{T}}$. Коэффициенты $C_{P1}^h$ и $C_{P1}^C$ равны $0.46/{\text{Pr}}_T^2$ и $0.46/{\text{Sc}}_T^2$, соответственно.

Слагаемое с перекрестной диффузией ${\mathrm{\mu }}_T\partial \stackrel{~}{{f^{\text{'}\text{'}}}^2}/\partial x_k\partial {\mathrm{\omega }}_f/\partial x_k$, которое появляется при переходе от уравнений класса ${\sigma }_f-{\epsilon }_f$ к уравнениям класса ${\mathrm{\sigma }}_f-{\mathrm{\omega }}_f$, подлежит удалению из моделей для частоты пульсаций. Дело в том, что это слагаемое ведет к нефизичному росту турбулентного числа Прандтля в слоях смешения и струях, что противоречит экспериментальным данным для сдвиговых течений [82].

Химические источники $Q_h^{\mathrm{chem}}=2\overline{{h^{\text{'}\text{'}}}_s{\dot{s}}_h}$ и $Q_C^{\mathrm{chem}}=2\sum _{j=1}^{N_{sp}}\overline{{Y^{\text{'}\text{'}}}_j{\dot{s}}_j}$ представляют собой корреляции между флуктуациями скаляра и их источниками, где ${\dot{s}}_h$ – скорость тепловыделения, ${\dot{s}}_j$ – источник массы для j-го компонента, они определяются через скорости реакций ${\dot{\omega }}_k$ и разность стехиометрических коэффициентов в левой и правой частях уравнения прямой реакции: $\Delta {\nu }_{jr}={\nu }_{jr}^f-{\nu }_{jr}^b$:

${\dot{s}}_j=-W_j\sum _{r=1}^{N_{react}}\left(\Delta {\nu }_{jr}{\dot{\omega }}_r^f-\Delta {\nu }_{jr}{\dot{\omega }}_r^b\right),{\dot{s}}_h=-\sum _{j=1}^{N_{sp}}\Delta h_j^0{\dot{s}}_j$

Здесь $W_j$ – молярная масса j-го компонента. В прямых реакциях $\Delta {\nu }_{jr}>0$ для исходных веществ и $\Delta {\nu }_{jr}<0$ для продуктов.

Существует множество подходов к моделированию корреляций источника скаляра вида $Q_h^{\mathrm{chem}}=2\overline{{h^{\text{'}\text{'}}}_s{\dot{s}}_h}$ и $Q_C^{\mathrm{chem}}=2\sum _{j=1}^{N_{sp}}\overline{{Y^{\text{'}\text{'}}}_j{\dot{s}}_j}$. Простейший подход, описанный в [41, 84], – алгебраическая модель, выражающая корреляцию  $2\overline{f^{\text{'}\text{'}}{\dot{s}}_f}$через средние источники, $\overline{{\dot{s}}_f}$, и среднеквадратичное отклонение пульсаций скаляра, $\sqrt{{\sigma }_f}$. Был проведен [85] анализ предварительно перемешанного турбулентного пламени на основе функций плотности вероятности (ФПВ), предложены аппроксимации данных слагаемых, основанные на переменной прогресса реакции. В текущей статье используется модель, основанная на методе PaSR для первого канала TCI:

$\sum _{j=1}^{N_{sp}}\overline{{Y^{\text{'}\text{'}}}_j{\dot{s}}_j}={\gamma }^{*}\sum _{j=1}^{N_{sp}}\left(Y_j^{*}-\widetilde{Y_j}\right){\dot{s}}_j^{*}, \overline{{h^{\text{'}\text{'}}}_s{\dot{s}}_h}={\gamma }^{*}\left(h_s^{*}-\widetilde{h_s}\right){\dot{s}}_h^{*}$ (4.11)

Это замыкание должно в основном оставаться положительным. Действительно, если j-й компонент представляет собой реагент, расходующийся при горении, то наиболее вероятно, что в тонких структурах массовая доля реагента меньше, чем в окружающем пространстве: $Y_j^{*}<Y_j^0$ и, следовательно, $Y_j^{*}<{\tilde{Y}}_j$ (см. (3.1)). В то же время, мы можем ожидать, что ${\dot{s}}_j^{*}<0$. Следовательно, произведение $\left(Y_j^{*}-\widetilde{Y_j}\right){\dot{s}}_j^{*}$ будет в основном положительным. Легко проверить, что оно также будет преимущественно положительным, если j-й компонент представляет собой продукт, образующийся при горении. Следовательно, в большинстве случаев $Q_C^{\mathrm{chem}}=2\sum _{j=1}^{N_{sp}}\overline{{Y^{\text{'}\text{'}}}_j{\dot{s}}_j}>0$. Аналогичный вывод может быть сделан для $Q_h^{\mathrm{chem}}$, так как явная энтальпия смеси должна увеличиваться вследствие тепловыделения. Однако, отрицательное значение корреляции тоже возможно, если состав и явная энтальпия в окружающем газе формируются не только вследствие обмена с тонкими структурами, но также подвержены сильному влиянию диффузии от соседних ячеек расчетной области. Чтобы предотвратить возможные осцилляции, которые не могут и не должны разрешаться в рамках RANS подхода, следует ограничить значение замыкания для $Q^{\mathrm{chem}}$ нулем снизу.

Граничные условия для модели PrOm сравнительно простые. Квадрат пульсаций скаляра на стенке принимается равным нулю ${\left.{\mathrm{\sigma }}_h\right|}_W=0$, ${\left.{\mathrm{\sigma }}_C\right|}_W=0$. Градиенты турбулентных чисел Прандтля и Шмидта предполагаются равными нулю на стенке: ${\left.\partial {\text{Pr}}_{\text{T}}^h/\partial n\right|}_W=0$, ${\left.\partial {\text{Sc}}_T/\partial n\right|}_W=0$. Эти условия совместно с соотношением (4.5) позволяют связать граничные значения частоты пульсации скаляра с граничным условием для $\mathrm{\omega }$.

Было также обнаружено [86], что модель градиентной диффузии Буссинеска может нарушать критерий реализуемости турбулентного потока. Критерий реализуемости формулируется согласно неравенству Шварца:

$\stackrel{~}{f^{\text{'}\text{'}}{u^{\text{'}\text{'}}}_i}\stackrel{~}{f^{\text{'}\text{'}}{u^{\text{'}\text{'}}}_i}\le 2\tilde{k}{\mathrm{\sigma }}_f$

Принимая во внимание аппроксимацию Буссинеска $\stackrel{~}{f^{\text{'}\text{'}}{u^{\text{'}\text{'}}}_i}\approx -\frac{\tilde{k}}{{\mathrm{\omega Pr}}_T^f}\frac{\partial \tilde{f}}{\partial x_i}$ и соотношение (4.1), получаем нижнюю границу для ${\text{Pr}}_{\text{T}}^f$:

$P{r}_T^f\ge \frac{1}{\mathrm{\omega }}\sqrt{\frac{1}{2}\frac{\tilde{k}}{{\mathrm{\sigma }}_f}\frac{\partial \tilde{f}}{\partial x_j}\frac{\partial \tilde{f}}{\partial x_j}}$

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

БЛАГОДАРНОСТИ

Авторы благодарят рецензента за внимательное прочтение работы и сделанные замечания и предложения, учтенные при доработке статьи. Представленные в статье новые результаты получены при финансовой поддержке Министерства образования и науки РФ (Договор № 14.G39.31.0001 от 13 февраля 2017).

Приложение. Для вывода алгебраической формулы для турбулентного числа Прандтля необходимо рассмотреть уравнение для турбулентного потока скаляра $\stackrel{~}{f^{\text{'}\text{'}}{u^{\text{'}\text{'}}}_i}$. Незамкнутое уравнение записывается так:

$\begin{array}{l}\frac{\partial \left(\overline{\mathrm{\rho }}\stackrel{~}{{\mathrm{f}}^{\text{'}\text{'}}{{\mathrm{u}}^{\text{'}\text{'}}}_{\mathrm{i}}}\right)}{\partial \mathrm{t}}+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{j}}}\left[\overline{\mathrm{\rho }}\stackrel{~}{{\mathrm{f}}^{\text{'}\text{'}}{{\mathrm{u}}^{\text{'}\text{'}}}_{\mathrm{i}}}{\tilde{\mathrm{u}}}_{\mathrm{j}}+\overline{\mathrm{\rho }}\stackrel{~}{{{\mathrm{u}}^{\text{'}\text{'}}}_{\mathrm{j}}{{\mathrm{u}}^{\text{'}\text{'}}}_{\mathrm{i}}{\mathrm{f}}^{\text{'}\text{'}}}+\overline{{\mathrm{p}}^{\text{'}}\mathrm{f}}{\mathrm{\delta }}_{\mathrm{ij}}-\tilde{\mathrm{\mu }}\frac{\partial \stackrel{~}{{\mathrm{f}}^{\text{'}\text{'}}{{\mathrm{u}}^{\text{'}\text{'}}}_{\mathrm{i}}}}{\partial {\mathrm{x}}_{\mathrm{j}}}\right]=\\ \underset{\text{FI}}{\underbrace{-\overline{\mathrm{\rho }}\stackrel{~}{{\mathrm{f}}^{\text{'}\text{'}}{{\mathrm{u}}^{\text{'}\text{'}}}_{\mathrm{j}}}\frac{\partial \widetilde{{\mathrm{u}}_{\mathrm{i}}}}{\partial {\mathrm{x}}_{\mathrm{j}}}-\overline{\mathrm{\rho }}\stackrel{~}{{{\mathrm{u}}^{\text{'}\text{'}}}_{\mathrm{i}}{{\mathrm{u}}^{\text{'}\text{'}}}_{\mathrm{j}}}\frac{\partial \tilde{\mathrm{f}}}{\partial {\mathrm{x}}_{\mathrm{j}}}}}+\underset{\text{FII}}{\underbrace{\overline{{\mathrm{p}}^{\text{'}}\frac{\partial \mathrm{f}}{\partial {\mathrm{x}}_{\mathrm{i}}}}}}-\underset{\text{FIII}}{\underbrace{2\tilde{\mathrm{\mu }}\overline{\frac{\partial {\mathrm{f}}^{\text{'}\text{'}}}{\partial {\mathrm{x}}_{\mathrm{j}}}\frac{\partial {{\mathrm{u}}^{\text{'}\text{'}}}_{\mathrm{i}}}{\partial {\mathrm{x}}_{\mathrm{j}}}}}}+\underset{\text{FIV}}{\underbrace{\overline{{{\mathrm{u}}^{\text{'}}}_{\mathrm{i}}{\dot{\mathrm{s}}}_{\mathrm{f}}}}}\end{array}$ (П.1)

Первое слагаемое в правой части, обозначенное как $\text{FI}$, является производством градиентами среднего поля скорости и градиентом среднего поля скаляра, эти слагаемые не требуют замыкания. Слагаемые, отмеченные $\text{FII}$, являются корреляцией пульсации давления с градиентом скаляра, которое тоже обычно разделяют на медленную и быструю часть. Первая часть ведет себя как медленная изотропизация турбулентного потока, которая похожа на диссипацию и пропорциональна самому потоку. Вторая часть члена FII более сложная, она также считается быстрой частью обменного слагаемого и представляет собой реакцию направления турбулентного потока на пульсации давления. Физическая интерпретация этого слагаемого заключается в его стремлении повернуть вектор турбулентного потока скаляра в направлении, параллельном направлению собственного вектора бездивергентного тензора скоростей деформации $S_{ij}^{dil}={\tilde{S}}_{ij}-\frac{1}{3}\frac{\partial {\tilde{u}}_m}{\partial x_n}{\delta }_{ij},$ который соответствует собственному числу, приводящему к разрежению ${\alpha }^S$, которое положительно и удовлетворяет условию ${\alpha }^S>{\mathrm{\beta }}^S>{\gamma }^S$. Последнее собственное число ${\gamma }^S$ отрицательное, тогда как знак ${\mathrm{\beta }}^S$ определяет форму трехмерного сжатия-растяжения жидкой частицы. Подобные рассуждения, анализирующие топологию градиентов скорости, использовались для описания производства скалярной диссипации в [45]. Простое замыкание для члена $\text{FII}$ записывается так:

${\text{FII}}_i=-C_{\Phi 1}\overline{\rho }\frac{1}{T_f}\stackrel{~}{f^{\text{'}\text{'}}{u^{\text{'}\text{'}}}_i}+C_{\Phi 2}\overline{\rho }{S}_{ij}^{dil}\stackrel{~}{f^{\text{'}\text{'}}{u^{\text{'}\text{'}}}_j}$ (П.2)

Здесь $C_{\mathrm{\Phi }1}=3.0$ и $C_{\mathrm{\Phi }2}=1.0;$ $T_f$ – смешанный масштаб времени, $T_f^{-1}\approx \sqrt{{\mathrm{\omega }}_f\mathrm{\omega }}.$ Конечно, это замыкание нужно лишь для вывода алгебраической модели. Существует множество работ, где представлены очень изощренные модели (см., напр., [86–88]). Простое замыкание (П.2), которое предлагается в этой статье, было протестировано на DNS базе данных турбулентного течения Рэлея–Тейлора с плавучестью [89] и продемонстрировало коэффициент корреляции, близкий к 0.4–0.45, тогда как гораздо более сложные модели не превышали значения коэффициента корреляции, равного 0.6. Коэффициент корреляции определяется следующим образом:

$\frac{{\displaystyle {\sum }_{ijk}^{N_{\mathrm{cells}}}{\text{FII}}_m^{\mathrm{real}}{\text{FII}}_m^{\mathrm{model}}}}{\sqrt{{\displaystyle {\sum }_{ijk}^{N_{\mathrm{cells}}}{\left({\text{FII}}_m^{\mathrm{real}}\right)}^2{\sum }_{ijk}^{N_{\mathrm{cells}}}{\left({\text{FII}}_m^{\mathrm{model}}\right)}^2}}}$

Процедура калибровки, похожая на градиентный спуск, описана в статье [90]. Согласно полученным коэффициентам корреляции упрощение, вносимое уравнением (П.2), можно считать приемлемым.

Следующее слагаемое, обозначаемое как $\text{FIII}$, является анизотропной диссипацией. Оно может быть описано более точно с дополнительным слагаемым, которое является сверткой потока скаляра с тензорами напряжений Рейнольдса:

${\text{FIII}}_i=\overline{\mathrm{\rho }}{\mathrm{\omega }}_f\left(C_{\in 1}\stackrel{~}{f^{\text{'}\text{'}}{u^{\text{'}\text{'}}}_i}-C_{\in 2}\frac{\stackrel{~}{{u^{\text{'}\text{'}}}_i{u^{\text{'}\text{'}}}_j}}{\tilde{k}}\stackrel{~}{f^{\text{'}\text{'}}{u^{\text{'}\text{'}}}_j}\right)$ (П.3)

Значения коэффициентов $C_{\in 1}$ и $C_{\in 2}$ были подобраны для оптимального описания базы данных DNS турбулентного течения Рэлея–Тейлора [89] и смогли воспроизвести точное значение ${\text{FIII}}_i$ с коэффициентом корреляции 0.73. Эффекты мелкомасштабной анизотропии в этом конкретном потоке были связаны с эффектами плавучести, что привело к появлению второго дополнительного слагаемого. Значения $C_{\in 1}$ и $C_{\in 2}$ оказались низкими, поэтому членом ${\text{FIII}}_i$ можно пренебречь при создании алгебраической модели. Легко продемонстрировать, что этот член пренебрежимо мал для изотропного случая, поскольку для изотропной турбулентности корреляция между градиентами скаляра и градиентами скорости определяется мелкомасштабными вихрями, для которых не существует выделенного направления градиента скаляра, поэтому мелкомасштабная скорость и скалярные градиенты статистически независимы, что приводит к нулевому пределу этого члена при числе Рейнольдса, стремящемся к бесконечности. Более того, этот член на самом деле является вектором и не должен иметь выделенного направления для мелкомасштабной изотропной турбулентности, поэтому им пренебрегают в случае пассивного скаляра. Однако он может быть существенным из-за мелкомасштабной анизотропии, которую можно наблюдать для активного скаляра в факельном режиме горения, где число Карловица $\mathrm{Ka}={\mathrm{\tau }}_{\mathrm{chem}}/{\mathrm{\tau }}_K$ меньше единицы. Рассматривая режим расширенного фронта пламени, где $\mathrm{Ka}\gg 1$ и много микромасштабных вихрей присутствует внутри фронта пламени, можно допустить локальную изотропность и пренебречь членом ${\text{FIII}}_i$.

Слагаемое $\text{FV}$ не существует для пассивного скаляра. Его замыкание для активного скаляра в режиме предварительно перемешанных флеймлетов рассмотрено Либби в [85], где предполагается Гауссова функция плотности вероятности скорости и добавляя ее зависимость от переменной прогресса реакции. Замыкание для этого источника на основе подхода PaSR, аналогичное (4.11), может быть сформулировано следующим образом:

${\text{FIV}}_i=\overline{{u^{\text{'}}}_i{\dot{s}}_f}={\gamma }^{*}\left(u_i^{*}-{\tilde{u}}_i\right){\dot{s}}_f\left(T^{*},f^{*}\right)$ (П.4)

В рамках классического подхода PaSR, которое используется в настоящей статье, скорость газа в тонких структурах и в окружающем пространстве одинакова: $u_i^{*}=u_i^0={\tilde{u}}_i$. Тогда (П.4) дает ${\text{FIV}}_i=\overline{{u^{\text{'}}}_i{\dot{s}}_f}\approx 0$.

Однако для предварительно перемешанных пламен, используя некоторые идеи из статьи Либби [85], можно предложить более точную трактовку этого члена, основанную на PaSR:

${\text{FIV}}_i=\overline{{u^{\text{'}}}_i{\dot{s}}_f}=C_B{\gamma }^{*}\left(1-{\gamma }^{*}\right)\frac{{\delta }_f}{{\tau }^{*}}{n}_i{\dot{s}}_f\left(T^{*},f^{*}\right)$ (П.5)

Здесь $n_i$ – нормаль к фронту пламени, направленная в сторону горячей области, ${\delta }_f$ – толщина пламени, ${\mathrm{\tau }}^{*}$ – характерное время протекания газа через тонкие структуры. Это замыкание легко получить, предположив совместную ФПВ, которая зависит как от мгновенной скорости, так и от дискретной случайной величины $\widehat{\gamma }$, которая равна единице в тонких структурах и нулю в окружающем пространстве. Среднее значение скорости, которое присутствует в гауссовом распределении, зависит от $\widehat{\gamma }$ следующим образом:

${\tilde{u}}_i\left(\widehat{\gamma }\right)={\tilde{u}}_i^0+\mathrm{\theta }\left(\widehat{\gamma }-\frac{1}{2}\right)\frac{{\mathrm{\delta }}_f}{{\tau }^{*}}{n}_i, {\tilde{u}}_i^0={\tilde{u}}_i-{\gamma }^{*}\frac{{\text{δ}}_f}{{\tau }^{*}}{n}_i$

Здесь θ – функция Хэвисайда. Интегрируя ${u^{\text{'}}}_i{\dot{s}}_f$ с совместной функцией плотности вероятности $\text{PDF}\left(\overrightarrow{u},\widehat{\gamma }\right)$, можно получить аналитическое выражение (П.5). Процесс вывода аналогичен описанному в статье Либби [85].

Наконец, вернемся к выводу алгебраического выражения для турбулентного числа Прандтля. Будем рассматривать f как пассивный скаляр, тогда членами FIII и FIV можно пренебречь. Для получения решения необходимо приравнять сумму слагаемых FI и FII в правой части уравнения к нулю. Очевидно, что производство градиентами скорости полностью компенсируется быстрой частью обменного члена, если в обоих случаях брать бездивергентные тензоры скорости деформации. Таким образом, два оставшихся слагаемых – это производство градиентами скаляра и изотропизация:

$C_{\mathrm{\Phi }1}\overline{\mathrm{\rho }}{T}_f^{-1}\stackrel{~}{f^{\text{'}\text{'}}{u^{\text{'}\text{'}}}_i}=-\overline{\mathrm{\rho }}\stackrel{~}{{u^{\text{'}\text{'}}}_i{u^{\text{'}\text{'}}}_j}\frac{\partial \tilde{f}}{\partial x_j}$

Для случая изотропных пульсаций скорости ($\stackrel{~}{{u^{\text{'}\text{'}}}_1^2}=\stackrel{~}{{u^{\text{'}\text{'}}}_2^2}=\stackrel{~}{{u^{\text{'}\text{'}}}_3^2}=2\tilde{k}/3$) это выражение превращается в формулу Буссинеска , которую можно приравнять к формуле градиентной диффузии с турбулентным числом Прандтля:

$\stackrel{~}{{\mathrm{f}}^{\text{'}\text{'}}{{\mathrm{u}}^{\text{'}\text{'}}}_{\mathrm{i}}}=-\overline{\mathrm{\rho }}\frac{2}{{\mathrm{C}}_{\mathrm{\Phi }1}}\tilde{\mathrm{k}}{\mathrm{T}}_{\mathrm{f}}\frac{\partial \tilde{\mathrm{f}}}{\partial {\mathrm{x}}_{\mathrm{i}}}=-\frac{{\mathrm{\mu }}_{\mathrm{T}}}{{\text{Pr}}_{\mathrm{T}}^{\mathrm{f}}}\frac{\partial \tilde{\mathrm{f}}}{\partial {\mathrm{x}}_{\mathrm{i}}};{\mathrm{\mu }}_{\mathrm{T}}=\overline{\mathrm{\rho }}\frac{\tilde{\mathrm{k}}}{\mathrm{\omega }}$

Тогда турбулентное число Прандтля определяется следующим образом:

${\mathrm{Pr}}_{\mathrm{T}}^{\mathrm{f}}=\frac{{\mathrm{C}}_{\mathrm{\Phi }1}}{2}\frac{1}{{\mathrm{\omega T}}_{\mathrm{f}}}={\mathrm{C}}_{\mathrm{\Phi }}\sqrt{\frac{{\mathrm{\omega }}_{\mathrm{f}}}{\mathrm{\omega }}}$

作者简介

V. Vlasenko

TsAGI; MIPT

编辑信件的主要联系方式.
Email: vlasenko.vv@yandex.ru
俄罗斯联邦, Zhukovsky; Dolgoprudny

R. Balabanov

TsAGI; MIPT

Email: vlasenko.vv@yandex.ru
俄罗斯联邦, Zhukovsky; Dolgoprudny

Wencha Liu

MIPT

Email: vlasenko.vv@yandex.ru
俄罗斯联邦, Dolgoprudny

S. Molev

TsAGI

Email: vlasenko.vv@yandex.ru
俄罗斯联邦, Zhukovsky

V. Sabelnikov

TsAGI

Email: vlasenko.vv@yandex.ru
俄罗斯联邦, Zhukovsky

参考

补充文件