Mathematical features of numerical simulation of non-stationary flow around the model in the shock tube

全文:

1. Введение

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

В настоящее время активно разрабатывается концепция «Электронная Аэродинамическая Труба» [1], которая подразумевает сопровождение эксперимента с целью последующей коррекции и пересчета экспериментальных данных. Созданы «Цифровые двойники» нескольких промышленно важных АДТ и имеется обширный опыт их практического применения [2–4]. Задача решена в рамках системы уравнений RANS [5], замкнутой дифференциальной моделью турбулентности [6, 7] и дополненной специальными граничными условиями [3], моделирующими работу реальной АДТ, имеющей камеру давления, перфорацию, поддерживающее устройство, камеру смешения и т.д. В данной работе ставится задача подготовки математического базиса для разработки «Цифрового двойника» ударной АДТ [8], в которой реализуется существенно нестационарный процесс, определяемый скачком уплотнения, проходящим по ее тракту.

Нестационарная задача для системы уравнений Навье–Стокса, осредненной по Рейнольдсу (URANS), выписана в нескольких монографиях, например, [11, 12]. Легко заметить, что системы RANS и URANS в данной статье идентичны, но это видимое совпадение. В подходе URANS время имеет четкий физический смысл. В уравнениях RANS, время добавлено в виде математического параметра для организации процедуры установления решения. Существуют и другие процедуры установления, например, итерационные, которые не требуют времени. Отличия RANS и URANS в данной реализации заключаются также в граничных и начальных условиях, и, кроме того, в особенностях применения используемых численных методов. Как уже упоминалось выше, в данной реализации RANS происходит установление решения по времени. При этом промежуточное состояние при условии, что процесс сходится, значения не имеет, важен результат. Для уравнений URANS исследуется развитие решения во времени и каждый промежуточный временной слой является определяющим. Подходы RANS и URANS достаточно часто подвергаются критике. Это связано с тем, что при выводе этих уравнений применяется процедура осреднения по всем масштабам турбулентных движений. Нестационарные процессы в URANS описываются корректно только в случаях, когда их временные масштабы намного больше масштабов турбулентности. Это ограничивает частоты исследуемых явлений. Следует отметить, что, наряду со «стационарными» моделями турбулентности, существуют и нестационарные, например, [13]. По этой причине вопрос о выборе модели решается в каждом отдельном случае, исходя из постановки задачи и результатов тестирования.

Задачи, связанные с распространением и взаимодействием нестационарных ударных волн, вызывают определенные трудности. Их аналитическое решение возможно только в приближенной или линейной постановке [14]. Что же касается нелинейных уравнений, то в настоящее время наиболее перспективными являются численные подходы, основанные на применении ЭВМ. Наиболее удобными являются методы типа Годунова [15–17], в которых естественным образом объединены возможности моделирования подвижных разрывов. Рассматривают два подхода к организации расчета: c выделением [18, 19] и размазыванием нестационарных разрывов [20]. Случай с размазыванием называется «сквозным расчетом», который характерен тем, что вычисления ведутся без специальной обработки зон с разрывами, которые появляются автоматически в виде областей с сильными изменениями параметров потока. Подход полностью оправдывает себя, т.к. позволяет моделировать сложные двумерные и трехмерные течения. В настоящее время получены интересные результаты, например, по дифракции [21] ударной волны на клине и др.

В методах сквозного счета ударная волна размывается за счет схемной вязкости, которая «работает» в ту же сторону, что и физическая вязкость. Тем не менее, численная и физическая вязкости принципиально различны по своей природе. Так, физическая вязкость размывает скачки уплотнения естественным образом, при этом процесс размытия описывается уравнениями Навье–Стокса, а численная вязкость увеличивает размеры области размытия за счет «паразитной» диссипации. С измельчением расчетной сетки влияние численной вязкости уменьшается пропорционально порядку численной схемы, и решение стремится к физически обоснованному. В случае невязкой задачи в постановке уравнений Эйлера ситуация усложняется, т.к. физически обоснованным решением для скачка уплотнения является разрыв, который методом сквозного счета не воспроизводится. По этой причине измельчение расчетной сетки вместо улучшения решения может приводить к его ухудшению за счет появления искажений, обусловленных локальной неустойчивостью, например, «карбункул-эффект» [22]. С другой стороны, «теряется» понятие аппроксимации скачка уплотнения. Так, в работах [23, 24] показано, что схема повышенного порядка аппроксимации имеет лишь первый порядок сходимости за фронтом ударной волны. В работах [25, 26] предлагаются различные способы преодоления указанных проблем, однако их обобщение на сложные задачи оказывается затруднительным. При этом определяющее значение имеет вопрос качества численной схемы. По соображениям, изложенным в работах [27, 28], приоритет отдается монотонным (или почти монотонным) схемам повышенного порядка точности, основанным на методе Годунова. В схемах Годунова повышение порядка аппроксимации по пространству обычно достигается процедурой реконструкции данных. Простейшими являются кусочно-линейные реконструкции, которые позволяют достичь второго порядка точности. Ключевую роль в создании монотонной схемы на основе кусочно-линейной реконструкции сыграла работа [29], в которой была предложена функция выбора минимального по модулю значения градиента. Указанная функция позволила записать численную схему повышенного порядка точности и обойти «запрет» знаменитой теоремы Годунова [30]. В дальнейшем была предложена схема типа MC [31] (monotonized central-difference), которая обладает меньшей диссипацией по сравнению с [29], но, строго говоря, не является монотонной, а удовлетворяет условию TVD (Total Variation Diminishing). Увеличение точности расчета продолжилось по пути построения численных методов на базе реконструкций пятого порядка. В рамках данной статьи использованы только два подхода. Это WENO5 (Weighted Essentially Nonoscillatory) и MP5, детальное описание которых можно найти в работe [32].

При решении нестационарной задачи важную роль играет способ интегрирования системы уравнений по времени. Из простейших можно указать явный метод Эйлера [33] первого порядка точности. Его дополняет неявный метод с дуальным шагом [34] по времени, который позволяет увеличить скорость расчета и проходить через заранее оговоренные временные интервалы. Определенный интерес представляет двухшаговый метод второго порядка, который является классическим методом Хьюна [35]. Указанный метод базируется на схеме RK (Runge–Kuta) и обладает повышенной устойчивостью при использовании TVD схем. Методы этого типа также идентифицируются как SSP (Strong Stability Preserving) [36]. Так, двухшаговый метод второго порядка точности с ограничением на шаг по времени, определяемым числом CFL = 1 (Courant–Friedrihs–Lewy) обозначается аббревиатурой SSP22. Ограничение на CFL является необходимым условием устойчивости численного решения дифференциальных уравнений в частных производных. Иногда для краткости число CFL называется числом Куранта. Пятишаговый метод четвертого порядка точности в этой классификации называется SSP54. В соответствии с работой [37] четырехшаговый метод SSP четвертого порядка построить невозможно. По всей видимости вариант, предложенный в работе [38], является оптимальным. Ограничение на шаг по времени в этом случае определяется числом CFL = 2.5.

Для настройки численного метода и оценки его точности применительно к проблеме создания цифрового двойника ударной АДТ используются различные тестовые задачи. Они делятся на две группы, которые отличаются типом продвижения ударной волны по тракту трубы. К первой группе относятся тесты, моделирующие продвижение волны в невозмущенном потоке [27]. При этом возможно появление препятствий, например, уступов [39]. Во второй группе тестов ударная волна движется в возмущенном поле. Здесь следует отметить задачу о нестационарном взаимодействии нескольких ударных волн [40]. Для тестирования метода расчета в случае неравномерного поля хорошо подходит задача Shu–Osher [35]. Она содержит синусоидальное возмущение плотности перед «правым» скачком уплотнения задачи Римана, которое приводит к сильному изменению решения указанной задачи.

При тестировании особое внимание следует уделять начальным ошибкам, возникающим у границы расчетной области вследствие искусственного «размазывания» ударной волны. Этот тип ошибки известен и описан в литературе как «start-up errors» [31]. На границе вследствие распада произвольного разрыва [30] возникают три возмущения со скоростями распространения «и–с», «и» и «и+с», где и – скорость потока, а с – скорость звука. Из-за наличия схемной вязкости указанные возмущения «размываются». При этом зона искусственного «размытия» накладывается на «размытие» вследствие действия физической вязкости, как молекулярной, так и турбулентной. Это приводит к образованию объединенного фронта возмущения, из которого трудно выделить физическую составляющую. Для этого требуется проводить специальные исследования.

В данной статье описывается математическая постановка задачи и численная схема для ее решения. Обсуждаются результаты тестовых расчетов и делаются оценки возникающих погрешностей. Исследуется обтекание изолированного конуса в неравномерном потоке. Рассматривается обтекание этого конуса в условиях запуска ударной АДТ. Дается описание экспериментальной установки. Проводится сопоставление расчетных и экспериментальных данных и делается заключение.

2. Система уравнений и математическая постановка задачи

Система нестационарных уравнений Навье–Стокса, осредненных по Рейнольдсу (URANS), записывается в виде:

$\frac{\partial \mathbf{u}}{\partial t}+\frac{\partial {\mathbf{f}}_j^{Eu}}{\partial x_j}+\frac{\partial {\mathbf{f}}_j^{NS}}{\partial x_j}=\mathbf{s}$

$\mathbf{u}=\left[\begin{array}{c}\mathrm{\rho }\\ {\mathrm{\rho u}}_{\mathrm{i}}\\ \mathrm{\rho E}\\ {\mathrm{\rho p}}_{\mathrm{m}}^{\mathrm{t}}\end{array}\right],{\mathbf{\text{f}}}_{\mathrm{j}}^{\mathrm{Eu}}=\left[\begin{array}{c}{\mathrm{\rho u}}_{\mathrm{i}}\\ {\mathrm{\rho u}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{j}}+{\mathrm{p\delta }}_{\mathrm{ij}}\\ (\mathrm{\rho E}+\mathrm{p}){\mathrm{u}}_{\mathrm{j}}\\ {\mathrm{\rho p}}_{\mathrm{m}}^{\mathrm{t}}{\mathrm{u}}_{\mathrm{j}}\end{array}\right],{\mathbf{\text{f}}}_{\mathrm{j}}^{\mathrm{NS}}=\left[\begin{array}{c}0\\ {\mathrm{\tau }}_{\mathrm{ij}}+{\mathrm{\tau }}_{\mathrm{ij}\mathrm{t}}\\ ({\mathrm{\tau }}_{\mathrm{ij}}+{\mathrm{\tau }}_{\mathrm{t}}){\mathrm{u}}_{\mathrm{j}}+({\mathrm{q}}_{\mathrm{j}}+{\mathrm{q}}_{\mathrm{j}\mathrm{t}})\\ {\mathrm{J}}_{\mathrm{j}}\left({\mathrm{p}}_{\mathrm{m}}^{\mathrm{t}}\right)\end{array}\right],\mathbf{\text{s}}=\left[\begin{array}{c}0\\ 0\\ 0\\ \mathrm{S}\left({\mathrm{p}}_{\mathrm{m}}^{\mathrm{t}}\right)\end{array}\right]$

Здесь u – вектор консервативных переменных, ${\mathbf{f}}_j^{Eu}$ и ${\mathbf{f}}_j^{\mathit{\text{NS}}}$ – вектора потоковых переменных, s – вектор источниковых членов, $u_i$ – компоненты вектора скорости, ρ – плотность, p – давление, E – полная энергия единицы массы, ${\tau }_{ij}=-\mu \left[\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}\frac{\partial u_i}{\partial x_i}{\delta }_{ij}\right]$ – компоненты тензора вязких напряжений, ${\tau }_{ijt}=\frac{2}{3}k{\delta }_{ij}-{\mathrm{\mu }}_t\left[\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}\frac{\partial u_i}{\partial x_i}{\delta }_{ij}\right]$ – компоненты тензора турбулентных напряжений, k – кинетическая энергия турбулентных пульсаций, µ и µ₁ – коэффициенты молекулярной и турбулентной вязкости. $q_j$, $q_t$ – молекулярный и турбулентный поток тепла. Для замыкания системы используются уравнения состояния $p=\mathrm{\rho }RT$ и одна из двух моделей турбулентности SA [6] или SST [7].

Решается краевая задача с известными начальными и граничными условиями. Решение получается методом интегрирования по времени. В начальный момент во всем пространстве задаются параметры набегающего потока, которые могут быть функцией времени. На твердых, в общем случае, подвижных границах выполняется условие прилипания. На входе и выходе из расчетной области задаются параметры, определяемые внешними условиями и инвариантами Римана. Интегрирование по времени осуществляется с шагами, удовлетворяющими ограничениям по числу Куранта [30] (CFL < 1) в ядре потока. В ряде случаев возможно получение установившегося и не зависящего от времени решения. В тех случаях, когда установившееся решение вследствие физических факторов не получается, а выходит на периодический цикл, может проводиться осреднение.

3. Численная схема и практическая реализация задачи

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

Для лучшего соответствия геометрии аэродинамической модели и АДТ расчетная область делится на блоки, которые стыкуются друг с другом. Это позволяет аппроксимировать сложные геометрии и получать качественные численные решения. Схема исследуемой АДТ приведена на рис. 1. Она взята непосредственно из работы [8]. Пример разбивки геометрии на блоки, которые привязаны к элементам АДТ, представлен на рис. 2. Хорошо видны поверхности конуса и элементов аэродинамической трубы. При построении расчетной сетки учитывается, что все получаемые решения являются сеточно-зависимыми. По этой причине максимально подробно учитываются местные особенности течения, в частности, положение скачков уплотнения.

 

Рис. 1. Схема ударной аэродинамической трубы УТ-1М [8]

1 – тепловая камера, 2 – электрический подогреватель; 3 – диафрагменный отсек; 4 – сверхзвуковое сопло; 5 – рабочая часть; 6 – оптическое окно; 7 – вакуумная емкость

 

Рис. 2. Математическая модель аэродинамической трубы УТ-1М

1 – камера высокого давления, 2 – форкамера, 3 – сопло, 4 – рабочая часть, 5 – конус

 

Для построения расчетной схемы используется система координат, связанная с сеточными линиями. Предполагается, что имеется ненулевой определитель матрицы Якоби J. Переход между декартовой и сеточной системами координат осуществляется по формуле: ${\mathbf{V}}_{\mathrm{\xi }}=\mathbf{J}\cdot {\mathbf{V}}_{\mathrm{x}}$. После перехода все действия производятся вдоль сеточных линий по одномерным соотношениям. На рис. 3 приведен пример пятиточечного сеточного шаблона. Границы ячеек изображены пунктирными линиями. Значения функций в центрах ячеек – кружками. Центрам ячеек приписываются целочисленные индексы, граням (пунктиры) – половинные.

 

Рис. 3. Пример пятиточечного сеточного шаблона

 

Вводятся следующие обозначения:

$\begin{array}{l}u=u_j,u_l=u_{j-1},u_r=u_{j+1}\\ \Delta u_l=u-u_l,\Delta u_r=u_r-u\\ h=h_j,h_l=(h_j+h_{j-1})/\mathrm{2,}{h}_r=(h_j+h_{j+1})/2\end{array}$

Линейная интерполяция возвращает значение $u_{inter}=(u_l{h}_r+u_r{h}_l)/(h_l+h_r).$ Для вычисления производных используется тот же шаблон. В результате получается:

$\begin{array}{c}{u^{\text{'}}}_l=\Delta u_l/h_l,{u^{\text{'}}}_r=\Delta u_r/h_r\\ u^{\text{'}}=({u^{\text{'}}}_l{h}_r+{u^{\text{'}}}_r{h}_l)=(\Delta u_l{h}_r^2+\Delta u_r{h}_l^2)/(h_l{h}_r^2+h_r{h}_l^2)\\ \Delta u=h{u}^{\text{'}}\end{array}$

Вторая производная записывается как:

$u^{\text{'}\text{'}}=({u^{\text{'}}}_r-{u^{\text{'}}}_l)/h=(\Delta u_r{h}_l-\Delta u_l{h}_r)/\left(h_l{h}_r{h}_l\right)$, а кривизна – $\Delta u^{\text{'}}=h{u}^{\text{'}}$

Для конструирования улучшенных схем второго порядка аппроксимации применяются математические функции:

$\mathrm{minmod}(\mathrm{x},\mathrm{y})=\frac{1}{2}\left(\mathrm{sign}\right(\mathrm{x})+\mathrm{sign}(\mathrm{y}\left)\right)\min \left(\right|\mathrm{x}|,|\mathrm{y}\left|\right)$

$\mathrm{median}(\mathrm{x},\mathrm{y},\mathrm{z})=\mathrm{x}+\mathrm{minmod}(\mathrm{y}-\mathrm{x},\mathrm{z}-\mathrm{x})$

Но существует и другой подход. Так, для конструирования TVD схемы MUSCL (Monotonic Upstream-centered Scheme for Conservation Laws) [28] значение производной вычисляется из условия, которое несколько отличается от функции .

$u^{\text{'}}=\left\{\begin{array}{ll}\mathrm{0,}& \mathrm{если}{u^{\text{'}}}_l{u^{\text{'}}}_r\le 0\\ \min \left(q\min \right({u^{\text{'}}}_l,{u^{\text{'}}}_r),0.5({u^{\text{'}}}_l+{u^{\text{'}}}_r\left)\right),& если{u^{\text{'}}}_l>0\mathrm{и}{u^{\text{'}}}_r>0\\ \max \left(q\max \right({u^{\text{'}}}_l,{u^{\text{'}}}_r),0.5({u^{\text{'}}}_l+{u^{\text{'}}}_r\left)\right),& \mathrm{если}{u^{\text{'}}}_l<0\mathrm{и}{u^{\text{'}}}_r<\mathrm{0,}\end{array}\right.$

где q = 1.25.

Более точная схема WENO5 основана на пространственной реконструкции пятого порядка [41]. Если решение гладкое и не используется сглаживание, то на равномерной сетке она дает пятый порядок аппроксимации по пространству. При наличии разрывов используется шаблон, дающий гладкое решение, а веса остальных шаблонов считаются близкими к нулю. Вариант монотонизации для схем высокого порядка и, в частности, для схемы WENO5 описан в [32]. Авторы [42] рекомендуют использовать сглаживание [41] и монотонизацию [32] одновременно. Процесс монотонизации разбивается на два этапа. Определяется ограничитель первого порядка на разрывах, а далее производится расширение диапазона на экстремумах. Схема WENO5 первоначально выписана для равномерной расчетной сетки. В случае существенного изменения размеров ячеек сетки (например, на границе блоков) порядок аппроксимации схемы падает до первого. Для решения этой проблемы в работе [43] предложен алгоритм учета неравномерности.

Для интегрирования по времени применяется многошаговая процедура перехода со слоя в момент времени n на следующий слой n +1. В работе исследованы три метода, имеющие разный порядок точности. Явный одношаговый метод Эйлера первого порядка точности описан в работе [30]. Двухшаговый явный метод SSP22 второго порядка точности подробно описан в работе [35]. Пятишаговый метод четвертого порядка точности SSP54. Коэффициенты пятишаговой схемы подобраны в работе [38] и обеспечивают устойчивость расчета.

4. Тестирование качества схемы расчета в нестационарной постановке

В данном разделе рассматривается тестовая задача движения плоской ударной волны по неоднородной среде в постановке Shu–Osher [35]. В начальный момент времени (t = 0) фронт ударной волны расположен в сечении x = 1. Далее ударная волна распространяется по газу, в котором плотность распределена по синусоидальному закону. Задача решается в приближении полной системы уравнений Эйлера. Начальный разрыв задается параметрами:

$\begin{array}{ll}({\mathrm{\rho }}_l,u_l,p_l)=\left(\mathrm{3.857143,}\mathrm{2.629369,}10.33333\right)& \mathrm{при}x<1\\ ({\mathrm{\rho }}_r,u_r,p_r)=(1+\epsilon \sin 5x,\mathrm{0,}1)& \mathrm{при}x\ge 1\end{array}$

Рассмотрен случай ε = 0.2. При этом реализуется полная конфигурация задачи Римана [30] (два скачка уплотнения с контактным разрывом между ними, левый скачок при данных условиях является слабым). Расчет ведется вплоть до момента времени t = 1.8. Задача Shu–Osher [35] не имеет точного решения. По этой причине «эталонное» решение получается численно и в дальнейшем используется для изучения качества различных схем. Рассмотрены две сетки с числом ячеек 16000 и 64000 ячеек. Используется схема WENO5 с параметром сглаживания ε = 10^-6. Для интегрирования по времени применены двухшаговая схема SSP22 и пятишаговая SSP54. В обоих случаях шаг по времени выбирался, исходя из условия CFL = 0.5. Полученные решения представлены на рис. 4. Толстой сплошной линией показано распределение плотности, полученное при решении задачи Римана в классической постановке, а тонкой сплошной, соответственно, решение, полученное для задачи Shu–Osher на подробной сетке (сетка 1) с применением схемы (WENO5, SSP54). Маркеры в виде крестов соответствуют решению на сетке в 16000 ячеек (сетка 2), которое получено по схеме (WENO5, SSP22). Анализ показывает, что оба решения с применением схемы WENO5 с погрешностью порядка $0.1\%$ совпадают друг с другом. По этой причине маркеры очень точно ложатся на тонкую сплошную линию. Это показывает, что для данного теста решение на подробной сетке (сетке 1) по схеме (WENO5, SSP54) с избытком заменяет отсутствующее точное решение. На графике можно выделить четыре области, соответствующие различным участкам течения. Область 1 располагается перед фронтом правого скачка уплотнения. Область 2 находится между фронтами правого скачка и контактного разрыва. Соответственно, область 3 занимает пространство между фронтами контактного разрыва и левого скачка уплотнения и, наконец, область 4 соответствует невозмущенному потоку перед левым скачком. На рис. 4 наблюдается расщепление решения, обусловленное различными условиями прохождения правого скачка уплотнения и контактного разрыва через область синусоидального возмущения. В области 2 проявляются нерегулярные колебания различной амплитуды, а в области 3 видны ассиметричные «пилообразные» осцилляции. В качестве эталона в дальнейшем используется решение (WENO5, SSP54) на сетке 1.

 

Рис. 4. Выбор эталонного решения для задачи Shu–Osher

 

Для дальнейшего анализа рассматриваются области 2 и 3, см. рис. 5. Для проведения практических расчетов выбирается сетка, которая содержит 400 ячеек (сетка 3). Применяются схемы четвертого (WENO5, SSP54) и второго порядка (MUSCL, SSP22) точности. Сопоставление показывает, что в области 3 на пике осцилляций (сверху) схема четвертого порядка соответствует эталону с погрешностью не хуже, чем 0.2% (график для увеличения масштаба построен в узком диапазоне координаты 5 > X > 3.5). Схема второго порядка дает на этих пиках погрешность 2%. При этом в донной области (снизу) точность расчета по схеме второго порядка еще ниже. В области 2 ситуация ухудшается, и схема второго порядка дает неприемлемо высокие ошибки. По этой причине, она в дальнейшем не рассматривается.

 

Рис. 5. Сопоставление решений четвертого и второго порядка точности на грубой сетке 3 при CFL = 0.5

 

На рис. 6 представлены результаты, полученные по схеме четвертого порядка с различными числами CFL (многошаговая процедура Рунге–Кутта позволяет работать с числами CFL , превышающими единицу [38]). Рассматривается только область 2 в узком диапазоне координаты 6.2 > X > 5.6, где уровень ошибок выше, чем в других местах. Сопоставление результатов расчета позволяет сделать вывод, что при значениях CFL = 0.5 и CFL = 1.25 полученные решения совпадают друг с другом (маркеры на графике перекрываются). При CFL = 2.5 диссипация становится заметной. Это учитывается при проведении практических расчетов.

 

Рис. 6. Сопоставление решений четвертого порядка точности при разных значениях CFL на грубой сетке 3 возмущения течения на границе расчетной области

 

В качестве следующего теста рассматривается осесимметричное обтекание острого конуса с углом полураствора θ = 5°. Геометрия расчетной области представлена на рис. 7.

 

Рис. 7. Геометрия расчетной области

 

Внешние границы расположены достаточно близко к поверхности конуса, что требует учета возможного отражения скачка уплотнения. На левой границе происходит втекание сверхзвукового потока. Это означает, что все возмущения входят через эту границу внутрь расчетной области и не оказывают обратного влияния. Технически данное граничное условие реализуется путем жесткого задания параметров внешнего потока в дополнительных ячейках, прилегающих слева к ячейкам на входной границе расчетной области. На правой границе реализуется граничное условие сверхзвукового вытекания. Это означает, что все возмущения выходят из расчетной области во внешний поток, а обратного влияния внешнего потока на расчетную область не существует. Технически данное граничное условие реализуется путем создания слоя дополнительных ячеек, прилегающих к выходной границе. В этих ячейках параметры потока задаются методом линейной экстраполяции соответствующих параметров из ячеек, прилегающих к выходной границе. На внешних границах расчетной области выполняется условие непротекания, что косвенно моделирует стенку ударной аэродинамической трубы. На поверхности конуса также ставится условие непротекания. В плоскости симметрии задачи также условие симметрии. Расчет проводится в два этапа. Сначала задается начальное поле с параметрами потока ${\mathrm{\rho }}_{t=0}=0.161317$ кг/м³, $u_{t=0}=567.17$ м/с, Па, $p_{t=0}=5930.32$, $T_{t=0}=128.07\text{\hspace{0.33em}K}$ и $M_{t=0}=2.5$ устанавливается стационарное поле потока вокруг конуса. Затем на левой границе вводится синусоидальное возмущение числа Маха по временному закону $\mathrm{\Delta }M=0.1\sin \left(10\mathrm{\pi }{u}_{t=0}t\right)$. Частота колебаний выбрана исходя из того, чтобы на длине конуса размещалось не менее пяти пиков осцилляций. Остальные параметры на левой границе вычисляются из условия постоянства полного давления и температуры:

$p_{total}=p_{t=0}/p\left(M_{t=0}\right)=p_{t=0}{\left(1+\frac{\gamma -1}{2}{M}_{t=0}^2\right)}^{\gamma /(\gamma -1)}$

$T_{total}=T_{t=0}/T\left(M_{t=0}\right)=T_{t=0}\left(1+\frac{\gamma -1}{2}{M}_{t=0}^2\right)$

$M=M_{t=0}+\Delta M,$ $p=p_{total}p\left(M\right)$

$T=T_{total}T\left(M\right),$ $u=M/\sqrt{\gamma RT}$

Расчет проводится в секторе размером в одну ячейку до тех пор, пока устанавливается квазистационарное периодическое решение. Это происходит через промежуток времени $\mathrm{\tau }=0.005$ с. Пример картины возмущенного обтекания конуса приведен на рис. 8, на котором в серой палитре представлено поле плотности. Хорошо видны возмущения, приходящие с левой границы расчетной области. Физическая картина течения на рис. 8 в некотором смысле повторяет течение в области 2 из задача Shu–Osher (см. рис. 4). Действительно, если перейти в систему координат, связанную со скачком уплотнения, то окажется, что синусоидально возмущенный поток набегает на указанный скачок. При этом, проявляются нерегулярные колебания плотности различной амплитуды, что и наблюдается в трехмерной реализации на рис. 8. При этом, поле плотности за скачком уплотнения возмущено достаточно сильно. Кроме того, скачок отражается от верхней стенки, но не приходит на конус.

 

Рис. 8. Поле плотности в момент времени t = 0.05 с

 

Решение получено на сетке 288×72×1 ячеек в постановке (WENO5, SSP54). На рис. 9 представлено сопоставление полученного решения с эталоном. В расчетной области выделены 3 точки. Первая точка располагается перед головным скачком уплотнения (x = 0.3, y = 0.2), вторая – за головным скачком уплотнения (x = 0.6, y = 0.2), а третья – за головным и отраженным скачками уплотнения (x = 0.9, y = 0.2). Наибольший интерес представляет точка 2, т.к. это «рабочая зона» исследуемого потока. Именно в этой точке проведено сопоставление с эталоном. Имеется соответствие с точностью порядка 0.1% в пиках. Наличие скачка уплотнения от конуса не повлияло на точность расчета, проведенного на достаточно редкой сетке. Это принципиально отличает данный тест от теста Shu–Oshera, где скачок уплотнения искажается синусоидальным полем перед его фронтом.

 

Рис. 9. Сопоставление с эталоном при наличии синусоидального возмущения на границе расчетной области

 

На рис. 10 приведены графики поведения плотности во времени для всех трех точек. Решение получено на сетке 288×72×1 ячеек в постановке (WENO5, SSP54). Наличие отраженного скачка уплотнения качественно не повлияло на результат. Следует отметить, что плотность газа увеличивается за скачками уплотнения. Кроме того, наблюдается заметный сдвиг фазы из за движения волн.

 

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

 

На следующем шаге приведен пример развития во времени решения в задаче вязкого обтекания острого конуса. Для этого построена подробная сетка у поверхности конуса, которая за пределами пограничного слоя плавно стыкуется с первоначальной «невязкой» сеткой. Размер первой ячейки равняется y⁺ = 1. Оценка размеров ячейки сделана на основании параметров пограничного слоя, полученных по эмпирическому методу Авдуевского [44]. Полный расчет вязкого течения выполнен в приближении URANS, описанном в разд. 1 данной статьи схемой WENO5 SSP54 с ограничением CFL = 1. Режим обтекания задан следующими параметрами потока на левой границе: слева от границы (в дополнительных ячейках) – M_Л = 2.5, P_{0 Л} = 810600 Па, T_{0 Л} = 293 K; справа от границы (в ячейках расчетного поля): M_П = 0, P_{0 П} = 101325 Па, T_{0 П} = 293 K. В результате взаимодействия начальных потоков образуются два скачка уплотнения, разделенные контактным разрывом. Указанная конфигурация движется направо в сторону установленного конуса, см. рис. 11,а.

 

Рис. 11. Стадии установления решения при сверхзвуковом обтекании конуса вязким потоком газа

а – скачки уплотнения и контактный разрыв в начальный момент времени;

б – появление отрыва в носовой части конуса;

в – развитый отрыв на поверхности конуса;

г – установившееся безотрывное обтекание конуса

 

Статическое давление на контактном разрыве (центральная стрелка на рис. 11,а) равняется P = 255124 Па. Перепад статического давления на правом скачке уплотнения равен $P_{\mathrm{К}.\mathrm{Р}.}/P_{\mathrm{П}}=2.52$, а на левом, соответственно, $P_{\mathrm{К}.\mathrm{Р}.}/P_{\mathrm{П}}=5.38$. После того, как левый скачок уплотнения в момент времени t = 0.003 с «набегает» на конус, появляется отрыв, см. рис. 11,б. Он имеет классическую «Лямбда» форму и начинается от носика. Контактный разрыв искривляется, т.к. поле потока за коническим скачком неоднородно. Правый скачок уплотнения уходит за пределы расчетной области. Отрывная зона растет и в момент времени t = 0.008 с занимает всю поверхность конуса, см. рис. 11,в. При этом, хорошо видна зона возвратного течения, точка присоединения потока и внутренние скачки уплотнения. В следующий момент времени зона отрыва уходит за пределы расчетной области и устанавливается сверхзвуковое обтекание конуса, см. рис. 11,г.

Для оценки времени установления решения в невязком ядре выбрана прямая линия с углом наклона Θ = 14°, что приблизительно соответствует среднему положению между поверхностью конуса и скачком уплотнения, см. рис. 11,г. На этой линии отмечены три точки. Известно, что с точностью до влияния пограничного слоя течение должно удовлетворять условию коничности. Другими словами, параметры потока в указанных точках должны совпадать друг с другом. На рис. 12 показано, как это условие выполняется во времени. Так, при t > 0.013 с относительные значения плотности с точностью лучше, чем 0.1%, совпадают друг с другом. Это соответствует моменту времени, когда отрывная зона уходит за пределы расчетной области.

 

Рис. 12. Выполнение условия коничности во времени

 

Несмотря на установление потока в невязком ядре, процесс установления в вязком потоке у стенки продолжается. На рис. 13 приведены профили пограничного слоя в три момента времени, построенные в сечении 3. Сопоставление показывает, что в момент времени t = 0.013 с профиль скорости u в пограничном слое остается «ненаполненным». Создается впечатление, что толщина пограничного слоя завышена. Но это впечатление ошибочно, так как имеет место неустановившееся течение. Процесс установления проходит достаточно быстро и при t > 0.023 с пограничный слой приобретает окончательную форму и перестает изменяться. С этого момента течение считается полностью установленным.

 

Рис. 13. Установление решения в пограничном слое

 

5. Моделирование нестационарного обтекания конуса в ударной АДТ

В данном разделе представлены результаты расчетно-экспериментального исследования. Эксперименты проведены в ударной аэродинамической трубе, работающей по схеме Людвига. Подробное описание указанной установки приведено в статье [45]. Исследовано обтекание конуса с углом 10° и длиной 0.96 м, см. рис. 14.

 

Рис. 14. Схема расположения приемников статического давления на поверхности конуса

 

Площадь дна конуса составляла 0.022 м², а технологический радиус носика – 0.08 мм. Конус устанавливался на шестикомпонентных внутримодельных тензометрических весах. Весы закреплялись в хвостовой державке, которая, в свою очередь, крепилась к вертикальной стойке. На поверхности модели имелись четыре приемных отверстия для измерения статического давления. Диапазоном измерения составлял ±35 кПа. Схема расположения приемных отверстий приведена на рис. 14. Три отверстия (Р1 – Р3) расположены на боковой поверхности конуса и одно (Р4) – на донной поверхности.

Наряду с экспериментальным исследованием, задача решена численно в упрощенной постановке. Конус размещен в ударной аэродинамической трубе. Математическая модель трубы достаточно точно соответствует оригиналу [8, 45]. Диаметр выходного сечения сопла составляет 0.5 м. Диаметр камеры высокого давления – 0.07 м, ее длина – примерно 12 м. Блок диафрагм находится на границе между «камерой высокого давления» и форкамерой. В начальный момент времени значение скорости во всем расчетном поле задано равным нулю. В «камере высокого давления» величина полного давления задана равной p₀ = 3203.1 кПа, а полная температура, соответственно, T₀ = 530.8 K В форкамере, сопле и рабочей части АДТ параметры потока заданы равными величинам p₀ = 150 Па, T₀ = 293.15 K.

Расчеты проведены по схеме (WENO5, SSP54) с условием CFL=1. Построена неструктурированная сетке, состоящая из 107721 ячеек, которая сгущена к твердым поверхностям с выполнением условия y⁺ = 1. В «камере высокого давления» сетка строится максимально грубой, т.к. поток в этой части экспериментальной установки не изучается. На верхней границе камеры ставится условие скольжения, хотя система уравнений содержит вязкие члены. Во всех остальных зонах АДТ на твердых поверхностях выполняется условие стенки с прилипанием и постоянной температурой T_W = 293.15 K. На «выходе» из расчетной области (справа) задается условие сверхзвукового вытекания потока.

Нестационарная картина обтекания конуса имеет достаточно сложную структуру. Область возмущенного течения, формирующего поток в рабочей части АДТ, достигает конуса в момент времени t = 0.0018 с и воздействует на него вплоть до момента t = 0.0055 с. При этом образуются скачки уплотнения и области завихрения. В качестве примера на рис. 15 приведены линии постоянства давления (изолинии) в момент времени t = 0.003 с.

 

Рис. 15. Физическая картина обтекания конуса в момент времени t = 0.003 с

 

Видно, что в указанный момент времени основной фронт ударной волны уже прошел лобовую поверхность конуса и взаимодействует с его дном. Скачки уплотнения, образующиеся в сопле АДТ, достигают носовой части конуса и пересекаются с отраженными скачками, а также с зарождающимся коническим скачком уплотнения. Описанная картина взаимодействия подтверждается показаниями датчиков статического давления. На рис. 16,а приведены показания датчика Р1. Значения статического давления отнесены к полному давлению в «камере высокого давления». Хорошо виден пик в момент времени t = 0.003 с, когда скачки уплотнения от сопла приходят на поверхность конуса в окрестности расположения датчика Р1. Далее следует провал с последующим восстановлением вплоть до времени t = 0.007 с. Расчет хорошо разрешает мелкие волны, которые видны на графике в виде осцилляций. В эксперименте указанные волны также видны, но «смазаны». Проведенная предварительно методическая работа позволяет доверять результатам расчета.

 

Рис. 16. Сопоставление расчетных и экспериментальных значений статического давления в измерительных точках 1, 3 и 4

 

Показания датчика Р3 качественно подтверждают описанную ранее картину течения, см. рис. 16,б. Отмеченный выше пик давления смещается во времени до t = 0.0035 с, что связано с его расположением на конусе. Наиболее сложная зависимость давления от времени соответствует датчику Р4 и представлена на рис. 16,в. Датчик расположен в донной области, где имеется отрыв потока и застойная зона. Используемый в данной работе метод на основе уравнений URANS не позволяет разрешать донные отрывные течения, характеризующиеся высокой степенью турбулентности потока. Поэтому к полученному результату в этой области следует относиться, как к модельному.

Заключение

Предложенная схема расчета на основе подхода (WENO5, SSP54) с условием CFL = 1 позволяет моделировать нестационарное обтекание конуса при нестационарных начальных и граничных условиях, включая запуск ударной АДТ.

Применение расчетных схем на основе подхода (MUSCL, SSP22) с условием CFL = 0.5 приводит к «смазыванию» физических осцилляций решения за счет повышенной диссипативности численной схемы и влечет за собой значительные погрешности. Она не рекомендуется к моделированию подобных задач.

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

Исследование выполнено при поддержке Российского научного фонда, грант № 23-19-00041.

作者简介

S. Bosnyakov

Central Aerohydrodynamic Institute named after prof. N.E. Zhukovsky

编辑信件的主要联系方式.
Email: bosnyakov@tsagi.ru
俄罗斯联邦, Zhukovsky

A. Gorbushin

Central Aerohydrodynamic Institute named after prof. N.E. Zhukovsky

Email: gorbushin@tsagi.ru
俄罗斯联邦, Zhukovsky

S. Matyash

Central Aerohydrodynamic Institute named after prof. N.E. Zhukovsky

Email: bosnyakov@tsagi.ru
俄罗斯联邦, Zhukovsky

S. Mikhailov

Central Aerohydrodynamic Institute named after prof. N.E. Zhukovsky

Email: mikhaylov@tsagi.ru
俄罗斯联邦, Zhukovsky

参考

补充文件