CC BY
50
Том XL III
УЧЕНЫЕ ЗАПИСКИ ЦАГИ 2012
№ 3
УДК 519.63
МУЛЬТИОПЕРАТОРНЫЕ СХЕМЫ ДО 18-го ПОРЯДКА ТОЧНОСТИ С ПРИЛОЖЕНИЯМИ К ЗАДАЧАМ НЕУСТОЙЧИВОСТИ И АКУСТИКИ СТРУЙ
М. В. ЛИПАВСКИЙ, А. Д. САВЕЛЬЕВ, А. И. ТОЛСТЫХ, Д. А. ШИРОБОКОВ
Описывается методика построения аппроксимаций и схем высокого порядка, основанная на использовании мультиоператоров, т. е. линейных комбинаций операторов из однопараметрических семейств. Рассматриваются варианты мультиоператоров, аппроксимирующих конвективные члены в уравнениях газовой динамики до 18-го порядка. Приводятся результаты численного моделирования неустойчивости дозвуковых и сверхзвуковых струй с генерацией акустических полей с помощью мультиоператорных схем и, в частности, численного моделирования скрич-эффекта в случае недорасширенных струй.
Ключевые слова: мультиоператоры, компактные схемы, уравнения Навье — Стокса, дозвуковые и сверхзвуковые струи, акустическое поле, неустойчивость струйных течений.
1. ВВЕДЕНИЕ
Одной из тенденций развития современных численных методов применительно к определенному классу задач, описываемых уравнениями математической физики, является повышение точности воспроизведения различных масштабов решений на больших интервалах времени. Ограничиваясь проблемами механики сплошных сред, к таким классам, имеющим практическое значение, можно отнести следующие задачи: моделирование атмосферных процессов;
применение прямого численного моделирования (DNS) и методики больших вихрей (LES) для описания турбулентности и других мелкомасштабных явлений;
моделирование генерации акустических полей, возбуждаемых струями реактивных двигателей, в интересах снижения уровня шума.
ЛИПАВСКИИ Михаил Витальевич
научный сотрудник ВЦ РАН
САВЕЛЬЕВ Александр Дмитриевич
кандидат физикоматематических наук, старший научный сотрудник ВЦ РАН
ТОЛСТЫХ Андрей Игоревич
доктор физикоматематических наук, профессор, заведующий отделом ВЦ РАН
ШИРОБОКОВ Дмитрий Алексеевич
кандидат физикоматематических наук, старший научный сотрудник ВЦ РАН
Традиционные способы повышения точности и разрешающей способности численных методов так или иначе связаны с увеличением числа базисных функций, на основе которых строятся аппроксимации исходных уравнений. В случае разностных схем это означает увеличение количества узлов в шаблонах (т. е. количества соседних узлов сетки, используемых в аппроксима-ционных формулах для каждого узла). Увеличение числа узлов, в свою очередь, порождает целый ряд неприятных проблем (в частности, возможность появления паразитических численных решений, необходимость специальных аппроксимаций в приграничных узлах, влияющих на устойчивость счета и т. д.). Хотя существуют способы борьбы с этими трудностями, при решении задач из перечисленных выше классов обычно точность схем не превышает 6-й порядок. Одним из примеров схем с многоточечными шаблонами, в котором улучшение свойства монотонности численных решений достигается путем использования линейных комбинаций аппроксимаций на различных шаблонах, являются схемы типа WENO [1].
Прямое численное моделирование процессов генерации звука является одной из сложных задач вычислительной аэродинамики и аэроакустики. Это связано не только с необходимостью выделять в решении акустические колебания давления с амплитудами, на несколько порядков меньшими уровней среднего давления, но и со способностью численных методов правильно описывать эволюцию коротковолновых составляющих численных решений в течение достаточно больших интервалов времени. Применительно к разностным методам это означает прежде всего использование схем высокого порядка, обеспечивающих малые фазовые и амплитудные ошибки в диапазоне длин волн, значительно превосходящих размеры шагов сеток. Однако высокий порядок, вообще говоря, не гарантирует малость этих ошибок в диапазоне коротких волн, поддерживаемых сетками, длина которых составляет несколько шагов сетки.
Чтобы сохранить в той или иной степени точность описания гармоник в широком диапазоне длин волн (или соответствующих им волновых чисел) требуются дополнительные средства, связанные с «настройкой» схем на хорошую разрешимость коротковолновых составляющих решений. Обычно это осуществляется в результате так называемой оптимизации схем путем управления свободными параметрами [2]. Эти параметры появляются как следствие расширения шаблонов без увеличения порядков аппроксимации.
Поясним приведенные выше соображения, сравнивая аппроксимации различных порядков пространственной производной в простейшем уравнении гиперболического типа
ди ди „
----+ a— = 0, a = const > 0
dt дх
с начальным условием в виде гармоники с волновым числом k
и( х, 0) = в1кх.
Такая задача обычно используется для оценки свойств разностных схем. Ее точным решением является волна, перемещающаяся слева направо с фазовой скоростью a :
и( х, t) = elk (x-at).
Если же заменить производную по х какой-либо разностной формулой на сетке с постоянным шагом h, то решение разностной задачи будет иметь вид:
u( х, t) = e“ d (kh)telk (x- a*(kh)t).
Оно характеризуется «схемной» фазовой скоростью a* (kh) и коэффициентом затухания гармоник d (kh) (для устойчивости схемы необходима неотрицательность функции d (kh)) в диапазоне волновых чисел, поддерживаемых сеткой, т. е. в диапазоне 0 < kh <п. Очевидно, что близость разностного решения к точному характеризуется отличием от единицы отношения a* / a и отличием от нуля d(kh), т. е. близостью к нулю фазовых и амплитудных ошибок.
2-і
а*/а, (І
2/ V
у /\\
за
Рис. 1. Зависимость фазовых (кривые 1, 4) и амплитудных ошибок (кривые 2, 3) от безразмерного волнового числа а = kh.
1, 2 — схема 1-го порядка; 3, 4 — схема 5-го порядка
На рис. 1 приведены зависимости этих ошибок от безразмерного волнового числа а = Ш для двухточечной аппроксимации первого порядка и аппроксимации пятого порядка из [3]. Цифрами на рисунке обозначены: 1, 4 — «схемная» относительная фазовая скорость; 2, 3 — коэффициент затухания гармоник.
Как видно на рис. 1, при малых Ш (длины волн >> h) ошибки малы в обоих случаях, однако диапазон их малости значительно шире для аппроксимации пятого порядка. Таким образом, чтобы правильно описать поведение гармоники с фиксированным волновым числом k , нужно проводить расчеты с таким шагом h, чтобы произведение Ш попало в диапазон малости ошибок. Если при этом задаться одинаковой погрешностью, то для некоторых k это может означать, что в случае схемы 1-го порядка требуются существенно меньшие значения h . В случае трехмерных нестационарных задач газовой динамики условие устойчивости приводит к увеличению числа
арифметических операций в т4 при уменьшении шага в каждом пространственном направлении в т раз. Это делает схемы низкого порядка малопригодными для решения некоторых задач даже при использовании самых мощных компьютеров.
Заметим, что фазовые и амплитудные ошибки, как можно показать, при а = Ш ^ 0 убывают как 0(ап), где п либо равно порядку, либо отличается от него на единицу. Однако это не гарантирует их уменьшения с увеличением п при больших значениях а <п. Если же требуется хорошее разрешение коротких волн, допускаемых сетками, как это происходит в случае задач аэроакустики, то необходимы дополнительные меры для расширения диапазона малости фазовых и амплитудных ошибок.
Мультиоператорные аппроксимации [4, 5] позволяют эффективно выполнить требования, предъявляемые к разностным схемам высокой точности и высокой разрешающей способности.
Во-первых, они могут иметь формально любой разумно высокий порядок.
Во-вторых, желаемое повышение порядка не связано с увеличением числа узлов в шаблонах и с усложнением процедуры вычислений; при этом, если использовать параллельные вычислительные системы с малыми временами обмена данными, затраты машинного времени практически могут оказаться слабо зависящими от порядков аппроксимации.
В-третьих, они зависят от параметров, которыми можно распоряжаться для оптимизации, направленной на высокое разрешение коротковолновых составляющих акустических спектров.
Пример применения мультиоператорных схем в задаче прямого численного моделирования неустойчивости сдвиговых слоев с разрешением всех масштабов содержится в [6]. Ниже приводится основная идея методики и ее последние реализации в некоторых задачах моделирования источников шума в струйных течениях.
2. МУЛЬТИОПЕРАТОРНЫЕ АППРОКСИМАЦИИ И СХЕМЫ
2.1. ОБЩАЯ ИДЕЯ МЕТОДА
Пусть имеется однопараметрическое семейство операторов (5), аппроксимирующих с порядком т на сетке с постоянным шагом к некоторый линейный оператор Ь. Зафиксируем М значений параметра 5 , 5 = 5І, і = 1,2,...,М. Для узла сетки с номером ] запишем разложение в ряд Тейлора действия Ь(5І )[и]оператора Ьк на достаточно гладкую функцию и(х) (т. е. результата применения этого сеточного оператора к функции и(х), рассматриваемой в узлах сетки в виде:
Ьк (5і )[и] , = [Ьи] , + От (5 і )кт + От+і(5 і )к
т+1
/ \і т+М -2 . і т+М-К -От+М-2 (5 і )к + 0(к ).
(1)
Здесь квадратные скобки означают, что функция непрерывного аргумента рассматривается в узлах сетки. При этом значения высших производных в узле х]- включены в коэффициенты
разложения. Оператор Ь может быть любым оператором, аппроксимируемым различными формулами численного анализа (интегралом, производной, оператором сдвига и т. д.)
Равенства (1), выписанные для 5 = 5, умножим на пока неизвестные коэффициенты уг-, такие, что
м
I
г=1
& = 1, (2)
а затем сложим. Учитывая (2), получим:
м мм
'„т+1
1гит+1(5г)
г=1 г=1 г=1
1угЬк (5 )[и]] = [Ьи]] +1угат (5 ^ +1Угат+1(5г )к”
(3)
м
I
г =1
Е, \7 т+М-2 . г\( 1 т+М-К
Угат+М-2(5 )к + 0(к У
Записав М -1 равенств, обращающих в ноль коэффициенты при кт,кт+1,...,кт+м 2 :
м м м
I Угат (5г к = 0 , IУгат+1(5г ^ = 0 ,- IУгат+м-2 (5 )к^^ = 0 (4)
г =1 г =1 г =1
и добавив к ним равенство (2), получим линейную систему для определения коэффициентов. Предположим, что она имеет решение при любом М. Тогда линейная комбинация операторов в левой части (3), названная в [4] мультиоператором, аппроксимирует исходный оператор Ь с порядком м + т -1:
м
Ьм (51,52,..., 5м) = !уЛ (5) = Ь + 0(кт+м-1). (5)
г=1
Очевидно, что в случае разрешимости системы мультиоператор существует и единственен для любого набора м различных значений параметра. Операторы Ь (5) будем называть базисными операторами. Таким образом, можно говорить, что повышение порядка в мультиоператор-ной методике достигается не путем увеличения числа базисных функций, а путем увеличения числа базисных операторов. Такой подход обладает следующим достоинством.
Пусть имеется вычислительная система с м процессорами (например, м — ядерный компьютер). Тогда вычисление действий мультиоператоров на известную сеточную функцию выглядит следующим образом. Эта функция раздается каждому процессору (или ядру), после чего процессоры (или ядра) одновременно производят вычисления действий базисного оператора со «своим» значением параметра. Затем результаты собираются и суммируются с уже известными коэффициентами. Это открывает возможность построения аппроксимаций желаемого порядка точности без каких-либо усложнений базисных операторов и затрат машинного времени. Эта процедура схематически изображена на рис. 2.
Основным моментом при построении мультиоператоров является разрешимость системы для коэффициентов при любых м . Она отсутствует в случае обычных разностных формул, но может иметь место в случае однопараметрических семейств компактных аппроксимаций оператора Ь . Простейшая форма компактных аппроксимаций имеет вид:
Ьк (5) = Ак (5)-1 Вк (5), (6)
где Ак (5) и Вк (5) — обычные разностные операторы, определяемые суммированием по узлам шаблона. Оператор Ак (5) предполагается легко обратимым. Это условие хорошо выполняется,
и
Вычисление Вычисление Вычисление
м
Вычисление /,дД-У, ■5Л/)=2У/^й
<=1
Рис. 2. Параллельное вычисление действия мультиоператора с использованием многопроцессорного или многоядерного компьютера
если он записывается на трех- или двухточечных шаблонах, что и предполагается в данной работе. Можно ожидать, что для операторов типа (6) коэффициенты при возрастающих степенях к в (1) содержат возрастающие степени 5, что приводит к разрешимости системы (2), (4). Потенциальная возможность быть источником базисных операторов при любых 5 требует исследования в каждом конкретном случае. Форма (6) открывает следующую привлекательную возможность создания формул численного анализа сколь угодно высокого порядка аппроксимации. В качестве Вк (?) рассматривается «стандартный» оператор (например, формула трапеций), а в качестве Ак (5) — некоторая его коррекция, зависящая от параметра 5. Затем фиксируются м различных значений 5 и строится мультиоператор, позволяющий получить порядки аппроксимации, пропорциональные м.
2.2. МУЛЬТИОПЕРАТОРНЫЕ АППРОКСИМАЦИИ ПРОИЗВОДНЫХ
2.2.1. Основной мотивацией развития мультиоператорного метода было построение высокоточных и эффективных схем для уравнений механики жидкости и газа. В этом случае существенную роль играет аппроксимация конвективных членов. Первые мультиоператоры использовали базисные операторы из семейств несимметричных («противопотоковых») компактных аппроксимаций 3-го порядка, возникших, по-видимому, впервые как обобщение предложенной конструкции такого типа [7], а также ее развитие в виде аппроксимаций 5-го порядка [8] (см. также [9]). Они имели вид поправок типа (6) к оператору несимметричной аппроксимации 1-го порядка, который на сетке с постоянным шагом к можно записать в виде:
Д(^) = (Д0 -5Д2)/2к, где Д0 и Д2 суть операторы центральных первых и вторых разностей:
Д^^и j — и j+1 и j 1, Д^и ^ — и ^+1 2и j + и j 1.
При 5 — 1 и -1 оператор Д(5) переходит, соответственно, в левые и правые двухточечные разности. При этом многообразие аппроксимаций производных можно представить в виде:
ди и ди
— - А-1 (Д(5) + Вк (5)); — * Д(5) + А~-1Вк (5), (7)
дх дх
где Ак — трехточечный оператор, а Вк — либо трехточечный оператор, либо суперпозиция трехточечных и обратных трехточечных операторов. Подробное описание этих семейств и основанных на них мультиоператоров можно найти в [6, 10]. Пары этих операторов, соответствующие 5 > 0 и 5 < 0 отличаются только знаком их диссипативной составляющей и позволяют аппрок-
симировать производные «против потока» в устойчивых схемах. Обозначим через Ьк (5) оператор в правой части (7), аппроксимирующий производные с порядком к = 3.5 . Разложения действий этих операторов имеют вид (1) с коэффициентами:
®т = Рі(5 ^ ат+1 = P2(s),■■■,
где Рі (5) суть полиномы степени і либо от 5 , либо от 5_1. Такой вид коэффициентов позволяет последовательно определять суммы
мм м
Еу5 , Е^А2, - ,Еу. і=1 і=1 і=1
м м м
Еъ5-1, Е^5-2 Ет
і=1 і=1 і=1
м-1 15і или
-м +1 і і
и свести матрицу системы (2), (4) к обратимой матрице Вандермонда. Решение этой системы может быть выписано в аналитическом виде. Оно обеспечивает существование и единственность мультиоператора (5), в котором
Ьь ^) = Ьк ^), к = 3.5. (8)
При этом можно показать, что Ьм (•$!,...,) и Ьм (-,^1,...,_,?и) имеют одинаковую (кососимметрическую) составляющую, которая аппроксимирует производную, а их самосопряженные составляющие имеют противоположные знаки. Чтобы эти составляющие были диссипативными, требуется выделить те значения параметров, для которых они положительны. Эта задача легко решается, если считать, что все sj распределены заданным способом между их минимальными и максимальными значениями 5Ш;П и 5шах. Тогда требуется найти эти значения таким образом, чтобы действительная часть Фурье-образа мультиоператора была положительной. Заметим, что эта функция получается суммированием с коэффициентами yj действительных частей Фурье-образов базисных операторов, выписываемых в аналитическом виде. При вычислении действия базисных операторов (и, следовательно, мультиоператоров) необходимо выполнить трехточечные прогонки.
2.2.2. В последнее время был предложен [11, 12] новый принцип построения однопараметрических семейств компактных аппроксимаций, в которых используются двухточечные операторы с существенно более экономичным обращением. Этот принцип годится для различных формул численного анализа, однако основное внимание уделялось аппроксимации конвективных членов в уравнениях Навье — Стокса. Двухточечные операторы, зависящие от параметра (обозначим этот параметр через с ), образуют пару «левого» и «правого» операторов N, Иг в зависимости от того, используется узел х]_1 или узел х^+1 при записи действия этих операторов в узле х]-. Эти операторы имеют вид:
"{иI = (1 + с)ы] _ сы]_1, Нгы^I = (1 + с)иу- _ си+1.
Простейшие операторы, имеющие структуру, аналогичную (7):
°' = 2%(4° + 3с 4 2" Зс"7'42) • = гК4"_ 3с 4 2 + 3с"_142) • (9)
также будем относить к категории левых и правых операторов. Эти операторы аппроксимируют производную с третьим порядком:
°к[и]I =Ш + °(^3), к =',г . (10)
Чтобы найти или , достаточно вычислить = Щх42иу- или = "_142иу-, решив систему = 42иу- или = 42иу- . В виду структуры операторов "' и Nг это легко
осуществляется методом бегущего счета. В первом случае стартовыми значениями являются граничные условия на левом конце интервала, а во втором — на правом конце.
Можно показать, что при с >_1/2 пара (Д(с),Вг(с)) аналогична паре (Ьк(я),Ьк(_5)) и может быть использована при построении противопотоковых схем.
Разложения в ряды действий операторов Б1 и Бг содержат полиномы от с с последовательно возрастающими степенями, причем коэффициенты при четных степенях к отличаются лишь знаками. Для иллюстрации сказанного приведем первые члены этих разложений:
Di [u ], =
Dr [u ] j =
-+cI h\X4) -12 б
u_x-|± + c |hU4) -12 б
1 c c
------1------1------
3O б б
1
2 | Л h4U X5)
j
2I
3O б б
v у
h4U X5>
O(h5),
O(h5).
(11)
Зафиксировав Мзначений сх,С2,...,См и использовав разложение типа (11), получим систему того же типа, что и в случае операторов Ьк (5). Нетрудно убедиться в том, что системы оказываются одинаковыми для Б1 и Бг. Их решения определяют «левый» и «правый» мультиоператоры
M
M
DM,l ^l — cM ) = Z YiDi (ci) , DM,r ^l — cM ) = Z YiDr (ci).
i=1
i=1
В дальнейшем предполагается, что ci распределены по некоторому закону между значениями cmin >-1/2 и cmax . Поскольку число обусловленности систем с матрицами Вандермонда быстро возрастает с M, целесообразно использовать распределение в виде нулей полинома Чебышёва для интервала [cmin, cmax], минимизирующее это возрастание. Такое же распределение использовалось в случае операторов Lk (5).
2.2.3. Описанные выше мультиоператоры целесообразно использовать, когда число значений параметров M < 5 и порядки не превосходят 9. При больших М величины max | уг- | оказываются слишком большими (порядка 103), что в некоторых случаях может привести к заметному влиянию на величину ошибок округления. Чтобы преодолеть это ограничение для получения рекордно высокой точности, в последнее время использовался следующий прием.
В качестве базисных операторов используем полусуммы
D(c) = (Di(c) + Dr(c))/2, c = c,,...,c
M
Разложения типа (3) теперь будут содержать только нечетные степени h и полиномы четных степеней p2k (c), k = 1,2,... При этом образующаяся система вида (2), (4) не будет сводиться к системам с матрицами Вандермонда и их обусловленность, как показала практика, слабо зависит от М. При том же числе М в случае разрешимости системы (2), (4) (ее нужно исследовать
в каждом конкретном случае) получается погрешность вида O(h1(M+1)), что существенно выше порядка погрешности при прежних базисных операторах. Таким образом, при M = 4 получаем порядок 10, а при M = 8 — порядок 18. Такой выбор соответствует, например, многоядерным процессорам Intel Core 2 Dual Q6600 и i7-920, которые можно использовать для параллельных вычислений. При этом свободные параметры cmin, cmax позволяют оптимизировать фазовые ошибки в диапазоне самых коротковолновых гармоник, поддерживаемых сетками.
Построенные таким образом мультиоператоры не содержат диссипативный механизм. Чтобы ввести его в схему, можно построить другой мультиоператор, использовав базисные операторы
D(ct) = (Di (c,) - Dr (ct ))/2, i = 1,...,M .
Порядок его действия будет равным O(h2M-1) . Однако для него нужно выбирать только те значения ct, которые обеспечивают его положительность.
2.2.4. Во всех расчетах на основе мультиоператорных схем использовались схемы с расщеплением потоков. Проиллюстрируем конструкцию таких схем на примере одномерной системы нелинейных уравнений вида:
*+*<!!> = O, (12)
дt дx
где u и f суть P -компонентные векторы. Будем предполагать, что наборы значений параметров
являются допустимыми в указанном выше смысле. Обозначим через D + тот мультиоператор из пары мультиоператоров, основанных на описанных выше базисных операторах, который
положителен. Тогда другой оператор из этой пары будет отрицательным. Назовем его D-. Используя эти обозначения, запишем схему для (12) в безындексном виде:
rs 1
— + Du = O, Du = — [ D+ (f(u) + Cu) + D~ (f(u) - Cu)], (13)
дt 2
где C = const > O, а затем представим ее как
rs 1 1
— + -(D + + D~ )f(u) + - C (D + - D~ )u = O.
дt 2 2
Кососимметричный оператор во втором слагаемом аппроксимирует производную по X
функции f, а положительный самосопряженный оператор в третьем слагаемом, умноженный на положительную постоянную, определяет диссипацию. Таким образом, получаем схему с положительным оператором, что обеспечивает ее устойчивость в приближении замороженных коэффициентов. Она удобна тем, что не требует вычисления собственных значений матрицы дГ /дu , как это часто делается при построении нецентральных схем. Кроме того, постоянная С может быть использована для управления интенсивностью диссипации. В частности, ее нулевое значение «отключает» диссипативный механизм схемы. В случае же C > O он способствует подавлению мелкомасштабных паразитических колебаний при решении конвективных и конвективно-диффузионных задач.
Схема (13) содержит аппроксимацию нелинейного члена в (12) независимо от способа аппроксимации производной по t . Использование тех или иных способов полной дискретизации (12) требует дополнительных оценок устойчивости. В случае нестационарных задач целесообразно использовать методы Рунге — Кутты того или иного порядка точности, оказывающие стабилизирующее влияние на процесс вычислений. В некоторых случаях оказывается возможным отключить диссипативный механизм мультиоператоров, положив C = O . В случае стационарных задач имеют смысл двухслойные неявные схемы со сравнительно легко обращаемыми операторами на верхнем слое.
Схема (13) легко обобщается на случай многомерных задач. Для этого достаточно использовать рассмотренные операторы для аппроксимации производных в каждом из пространственных направлений.
2.3. ОЦЕНКИ ТОЧНОСТИ В СЛУЧАЕ МОДЕЛЬНЫХ ЗАДАЧ
Для иллюстрации приведем несколько примеров сравнения с точными решениями. В следующей таблице приведены погрешности решения периодической задачи для одномерного уравнения Бюргерса в зависимости от числа узлов сетки N.
N WENO-5 Ь5 Ь59 Ь10 Ь18
Е к Е к Е к Е к Е к
16 32 1.3Е-2 1.2Е-3 3.4 6.1Е-3 4.4Е-4 3.8 1.3Е-03 6.6Е-06 7.7 1.3Е-3 6.6Е-6 7.7 1.3Е-03 7.6Е-06 7.4
64 9.5Е-5 3.7 1.6Е-5 4.7 5.4Е-09 10.0 5.4Е-9 10.3 1.4Е-09 12.5
128 3.3Е-6 4.8 3.3Е-6 5.0 5.6Е-12 9.9 4.9Е-12 10.1 2.1Е-14 16.0
В таблице приняты следующие обозначения: WENO-5 — схема типа WENO 5-го порядка; Ь5 — компактная схема 5-го порядка; Ь59 — схема с мультиоператором 9-го порядка, основанном на базисных операторах 5-го порядка вида (8); Ь10 и Ь18 — схемы с мультиоператорами 10-го и 18-го порядков соответственно, основанными на базисных операторах (9); Е — отличие численного решения от точного в норме С; к — локальные порядки сходимости.
Из таблицы следует, что в случае мультиоператоров точность, близкая к машинной и превосходящая на 7 порядков точность в случае других схем, достигается при вполне реальном на современном уровне числе узлов сетки. Это создает перспективу описывать акустическое давление на 5 — 6 порядков меньшее, чем среднее давление в потоке.
Следующий пример иллюстрирует способность мультиоператорной схемы правильно описывать распространение коротковолновых компонент на больших интервалах времени. Известная тестовая задача Тама [13] состоит в следующем. Волновой пакет (рис. 3), центр которого в начальный момент находился при х = 0 , описывается функцией:
и(0,х) = [2 + со8(ах)]ехр(- 1п(2)(х/10)2), а = 2.3 .
Решается уравнение и( + их = 0 с такими начальными данными и шагом Дх = 1. Требуется сравнить решение с точным в момент времени ї = 800 .
Расчеты проводились по оптимизированной мультиоператорной схеме с аппроксимацией 10-го порядка производной по х. Ее фазовые ошибки, характеризуемые отличием отношения фазовых скоростей а* /а от единицы, приведены на рис. 4. Там же представлена аналогичная кривая для схемы с аппроксимацией 18-го порядка. Сравнивая рис. 4 и рис. 1, можно заметить, что эти схемы обладают существенно лучшими разрешающими способностями (почти что спектральным разрешением), обеспечивая малость ошибок вплоть до значений а = кк = 2.7 и а = кк = 2.85 для схемы 10-го и 18-го порядков соответственно. Видно, что численное решение (отмеченное маркерами) визуально совпадает с точным (сплошная кривая).
Рис. 3. Решение уравнения переноса для волнового Рис. 4. Зависимость фазовых ошибок от безразмерного
пакета: волнового числа а = кк:
------------точное решение; • • • — численное решение 1 — мультиоператорная схема 10-го порядка; 2 — мультиопера-
торная схема 18-го порядка
3. ПРИЛОЖЕНИЕ К НЕКОТОРЫМ ЗАДАЧАМ АЭРОАКУСТИКИ СТРУЙ
Одним из важных направлений аэроакустики является поиск путей снижения шума от струй реактивных двигателей. При решении этой проблемы существенную роль играет численное моделирование процесса генерации акустических волн сдвиговыми слоями на границах струй, истекающих из сопл. Сдвиговые слои являются потенциально неустойчивыми, причем эта неустойчивость, согласно теории, первоначально развитой для уравнения Гинзбурга — Ландау, может быть двух типов. В первом случае струя может быть усилителем попавших в нее возмущений, которые сносятся вниз по потоку (конвективная неустойчивость), во втором — струя работает как осциллятор и генератор излучаемых звуковых волн, при этом возмущения могут распространяться во все стороны, в том числе и вверх по потоку. Реализация режимов зависит от параметров течения, в частности, от толщины сдвигового слоя и температуры внутри струи. При этом режим течения в струе определяется главным образом параметрами струи и окружающего пространства и, практически, мало зависит от начальных и внешних возмущений. Другими словами, струя обладает своими резонансными частотами [14].
Ниже приводятся результаты применения вариантов мультиоператорных схем при решении акустических задач, для некоторых из них существует возможность сравнения с экспериментальными данными. Предварительно заметим следующее.
В акустических экспериментах для струй, которые рассматриваются в данной работе и к которым часто обращаются акустики, исследуются главным образом акустические поля, а полная информация о газодинамических полях и условиях их образования в приводимых публикациях отсутствует.
С другой стороны, для точного воспроизведения процессов излучения звука струями в условиях эксперимента существенной является возможность промоделировать все течение. Ввиду отсутствия необходимых для этого данных приходится создавать некоторые «синтетические» струи, которые приближенно соответствуют экспериментальным. Такая методика, реализуемая путем построения некоторых аналитических аппроксимаций полей течения с использованием небольшого количества имеющихся в распоряжении данных и известных газодинамических закономерностей, применялась, например, в [15] и [16]. Синтетические струи используются также для исследования процессов развития неустойчивости с генерацией звуковых полей без привязки к конкретным экспериментальным данным [19].
Чтобы осуществить сравнение расчетных и экспериментальных данных, в настоящей работе также пришлось восстанавливать поля течений в струях, близкие к экспериментальным, но используя численное моделирование. При этом точность такого восстановления зависит не столько от точности самих расчетов, сколько от полноты информации об условиях экспериментов. Соответственно, от этого же зависит и точность воспроизведения акустических полей, измеряемых в экспериментах. Таким образом, здесь можно говорить о сравнении скорее на качественном, чем на количественном уровне.
3.1. ГЕНЕРАЦИЯ ЗВУКА ВОЛНАМИ НЕУСТОЙЧИВОСТИ В СВЕРХЗВУКОВЫХ СТРУЯХ
В качестве первой задачи рассматривается численное моделирование акустических полей в рамках теории Тама о волнах неустойчивости, распространяющихся вниз по потоку в сверхзвуковых струях с тонкими сдвиговыми слоями.
Излучение звука при развитии неустойчивости происходит следующим образом. Сдвиговые слои на начальном участке сверхзвуковой струи характеризуются большими поперечными градиентами скорости, что приводит к неустойчивости и росту амплитуд гармоник. Например, монохроматическая волна от какого-то внешнего источника на срезе сопла увеличивает свою амплитуду при распространении вдоль оси. Вниз по течению толщина сдвиговых слоев увеличивается, а градиент скорости уменьшается, что ведет к затуханию волны. Возрастающая, а затем убывающая амплитуда приводит к широкополосному спектру волновых чисел. В коротковолновом диапазоне этого спектра возможны гармоники со сверхзвуковой фазовой скоростью относительно окружающей среды, они и являются источником звука.
Мы используем метод линеаризации [15] уравнений Эйлера около поля среднего течения. Заметим, что обобщение этой теории на неизэнтропические сверхзвуковые струи и доскональное ее тестирование представлено в [16]. Фактически мы изменили аналитический метод [15], применяя трехмерный численный код для моделирования акустических полей в осесимметричных
струях. Другое отличие от [15] состоит в вычислении среднего течения, соответствующего эксперименту, а не использовании аппроксимаций экспериментальных данных. Результаты расчетов сравниваются с аналитической теорией и экспериментальными данными [17].
3.1.1. Вычисление среднего течения
Поле струи было получено расчетным путем на основе осредненных уравнений Навье — Стокса и дифференциальной двухпараметрической модели турбулентности [18]. Диаметр выходного сечения сопла, числа Маха и Рейнольдса, нерасчетность истечения и размеры расчетной области соответствовали исходным данным эксперимента [17]. В то же время отсутствовали данные о форме сопла и турбулизации истечения воздушной струи. Поэтому в расчетах использовалось коническое сопло. Значение кинетической энергии турбулентности в критическом сечении подбиралось таким образом, чтобы получить распределения чисел Маха вниз по течению за соплом, близкие к экспериментальным. Выбранное таким образом, оно составило 5% от квадрата скорости звука. Расчетная методика использовала компактные аппроксимации 5-го порядка [9].
Распределения числа Маха и относительного давления вдоль оси струи представлены на рис. 5 соответственно сплошными и штрихпунктирными линиями. Там же кружочками приведены экспериментальные значения числа Маха из [17]. Отличия на начальном участке струи вызваны разной формой геометрии сопл: конической в расчете и специально профилированной в эксперименте. На выходе из конического сопла формируется система скачков уплотнения, более сильная, чем в эксперименте. На расстоянии примерно 12 диаметров от среза скачки угасают, и ниже по течению результаты расчета соответствуют данным эксперимента.
Рис. 5. Распределения числа М (------) и относительного давления (— • — • —)
вдоль оси струи; • • • — экспериментальные значения числа М из [17]
у г
у/г у/г \ \ ^
а) , б> А \ в)
\ 1 - А Чч. і \ \
•ч •
• •
•
• \т • • \ ✓' •
• г / і Г ,■% \ / #\
• 2 * » » і 0 й ^ о
Рис. 6. Распределения числа Маха в сечениях струи (-------), профили кинетической энергии турбулентности, отне-
сенные к квадрату звуковой скорости в критическом сечении сопла (— • — • —), экспериментальные значения числа М из [17] (• • •). Расстояния от кромки сопла в диаметрах составляют 1 (а), 5 (б), 10 (в) и 15 (г)
Распределения числа Маха в сечениях струи на расстояниях 1, 5, 10 и 15 диаметров от кромки сопла представлены сплошными линиями на рис. 6 в порядке слева направо. Качественно они согласуются с экспериментом (маркеры). Однако получить меньшие различия между расчетом и экспериментом оказалось затруднительно из-за отсутствия необходимых исходных данных.
На тех же графиках штрихпунктирными линиями приведены профили кинетической энергии турбулентности, отнесенные к квадрату звуковой скорости в критическом сечении сопла. Турбулентность генерируется в сдвиговом слое. Вблизи кромки сопла за счет резких градиентов скорости возникает сильный рост энергии турбулентности. Ниже по потоку, где скорость струи гасится при смешении ее с окружающим воздухом, происходит постепенное разглаживание профиля кинетической энергии и снижение ее значений.
3.1.2. Вычисление акустических полей
С использованием поля среднего течения были получены линеаризованные уравнения Эйлера. При этом отклонения нестационарных возмущений газодинамических параметров от параметров среднего течения рассматривались как флуктуации. В отличие от расчетов среднего поля использовалась декартова система координат с «диссипативными» демпфирующими зонами вблизи выходной и боковых границ. На входной границе задавались возмущения давления вида:
Р(г, ф, t) = / (г) 8т(ю I + пф),
где ф — азимутальный угол, или моделировалось возмущение подобное белому шуму. В расчетах использовались мультиоператорные схемы 7-го порядка точности. Представленные результаты были получены на сетке 61x81x81 узлов со сгущением вблизи струи.
Изолинии флуктуаций давления для п = 0 и 1 приведены на рис. 7, а, б. Число Струхаля БЬ = 0.2 соответствует частоте возбуждения на входной границе. На рис. 7, в показаны изолинии флуктуаций давления для возбуждения типа белый шум, п = 0 .
Распределение флуктуации пульсационного потока массы (ри) в сдвиговом слое сравнивается с результатами эксперимента [17] на рис. 8 (п = 0, БЬ = 0.2 слева, белый шум справа). Амплитуда давления на входе выбиралась так, чтобы максимальные значения флуктуаций в расчетах и в эксперименте совпадали.
Расчетные кривые приблизительно совпадают по форме с экспериментальными зависимостями, несколько отличаясь от них значениями пульсаций. В виду приближенного воспроизведения среднего поля экспериментальной струи, а также относительной малости пульсационных величин, измеряемых экспериментально, можно говорить (с учетом сделанного выше замечания о соответствии расчетных и экспериментальных данных) об их удовлетворительном согласии.
На рис. 9 представлена диаграмма направленности для БЬ = 0.2 (давление приводится в децибелах), обычно используемая для указания уровня акустического давления для каждого пространственного направления. На рис. 9 экспериментальные данные в виде маркеров располагаются достаточно близко к расчетной кривой, указывая на отличия до 3 дБ.
Рис. 8. Распределение флуктуации пульсационного потока массы в сдвиговом слое:
---------результаты расчетов; □ □ □ — эксперимент [17]
Таким образом, данная расчетная методика, основанная на модели излучения звука сверхзвуковыми струями в рамках теории волн неустойчивости, не показывает существенных различий вычисленных параметров акустических полей и полученных в условиях рассматриваемого эксперимента. Отличия в диаграмме направленности порядка 3 дБ, полученные в результате лишь приближенного воспроизведения среднего поля турбулентной струи, представляются вполне разумными.
Отметим, что в виду применения пространственного акустического кода данная вычислительная методика пригодна для аналогичных расчетов в случае пространственных сопл — достаточно лишь использовать расчетные данные для среднего трехмерного течения. Заметим, однако, что при данном подходе струя служит усилителем шума на срезе сопла (являясь при этом широкополосным фильтром), но она не является генератором собственного шума. В следующем разделе рассматриваются акустические поля, излучаемые самой струей в силу ее неустойчивости.
3.2. ГОРЯЧИЕ ДОЗВУКОВЫЕ СТРУИ
Для прецизионного моделирования различных режимов течения горячих струй использовались уравнения Навье — Стокса и мультиоператорные схемы 9-го и 10-го порядков с обращением, соответственно, трехточечных и двухточечных операторов. Они позволяли описывать акустические возмущения без линеаризации уравнений или введения добавок к среднему течению (как это делалось в работе [19]). Задача решалась в осесимметричной постановке в координатах х — вдоль течения струи, г — радиус поперек струи. Линия г = 0 соответствует оси симметрии струи в начальный момент времени. Расчетная область представляет собой прямоугольник [0, хтах ] X [0, Гтах ] . В начальный момент времени в струе задавались профили скорости, плотности и температуры с учетом соотношений Крокко — Буземана, не зависящие от продольной координаты х . На входе в область (х = 0) и на внешней границе г = гтах для предотвращения отражения возмущений внутрь расчетной области задавались «мягкие», характеристические условия. На оси струи (г = 0) ставились условия симметрии. В окрестности выхода из области к уравнениям добавлялись вспомогательные конвективные слагаемые и ставились условия сноса (экстраполяции из расчетной области в граничные узлы сетки).
Рис. 9. Диаграмма направленности для БИ = 0.2 (на осях координат отмечены уровни звукового давления в дБ; -----------------расчет; □ □ □ — эксперимент)
Рис. 10. Изолинии завихренности в меридиональной плоскости струи. Момент слияния вихрей
Рис. 11. Пространственно-временные диаграммы давления:
а — M = 0.1, £ = 0.55, 0 = 25; б — M = 0.1, £ = 0.65, 0 = 25; в — M = 0.1, £ = 0.4, 0 = 10; г — M = 0.1, £ = 0.55, 0 = 25; д — M = 0.4, £ = 0.55,
0 = 25; е — M = 0.8, £ = 0.55, 0 = 25
Потеря устойчивости выражалась в образовании и конвекции вихревых колец (со спариванием или без) в асинхронном или синхронном режимах. В последнем случае струя играла роль осциллятора. При фиксированном числе Маха режим течения определяется двумя параметрами: отношением температуры во внешнем пространстве к температуре в центре струи в начальном поперечном сечении £ = Тх / Т и отношением радиуса струи к толщине сдвигового слоя на границе струи (толщине вытеснения импульса) 0 = Я / а. Меридиональное сечение вихревого поля с событием спаривания колец в фиксированный момент времени изображено на рис. 10. В левой части рисунка хорошо видно образование очередного одиночного вихря. В центре происходит слияние двух одиночных вихрей в «двойной» вихрь. У правой границы мы видим единый вихрь, получившийся в процессе слияния двух одиночных. На рис. 11 изображены пространственновременные диаграммы течений, использованные в [19] при различных значениях £ и 0 и фиксированном числе М = 0.1 (рис. 11, а — в) и различных числах М (рис. 11, г — е). Здесь представлены значения давления вдоль линии г = 1, т. е. линии, проходящей по первоначальной границе струи. Темные полосы представляют низкие значения давления, это указывает на центры вихрей. Рис. 11, а соответствует значениям £ = 0.55 , 0 = 25 . При этих значениях реализуется синхронный режим с регулярным слиянием пар вихрей. Одиночные вихри во все моменты времени образуются на одном и том же расстоянии от входа в расчетную область (х « 2), и на расстоянии х«5.5 одиночные вихри попарно сливаются, образуя вихри больших размеров. Изменение параметра £ до 0.65 при том же значении 0 = 25 (рис. 11, б) приводит к возникновению асинхронного режима и нерегулярному слиянию пар вихрей. Как и в первом случае, происходит
а — М = 0.1; б — М = 0.4; в — М = 0.8
объединение пар вихрей, но отсутствует регулярность по времени, и слияния случаются на разном удалении от входа в расчетную область. При значениях S = 0.4, 0 = 10 (рис. 11, в) наблюдается синхронный режим, образуются одиночные вихри без слияния. На рис. 11, г — е представлены пространственно-временные диаграммы для значений S = 0.55, 0 = 25 : г соответствует М = 0.1, регулярное слияние вихрей; д — М = 0.4, нерегулярное слияние вихрей; е — М = 0.8, тоже нерегулярное слияние вихрей, но в отличие от предыдущего случая можно заметить вторичное спаривание «двойных» вихрей. На рис. 12 приведены акустические поля в расчетной области для этих же чисел Маха в децибелах. При М = 0.1 (рис. 12, а) струя излучает звук довольно протяженной областью. Четко не просматриваются какой-либо центр излучения и преимущественное направление. При М = 0.4 (рис. 12, б) обозначается направление излучения, примерно под углом в 40° от направления истечения струи. При М = 0.8 (рис. 12, в) четко прослеживается центр излучения, находящийся около оси струи в месте около первичного слияния вихрей. Хорошо видны два направления излучения звука.
3.3. МОДЕЛИРОВАНИЕ СКРИЧ-ЭФФЕКТА
При истечении из сопл сверхзвуковых недорасширенных струй наблюдается так называемый скрич-эффект («screech»). Он состоит в том, что образованные скачками уплотнения ячейки («бочки») в процессе их колебаний генерируют акустические волны, распространяющиеся против потока. Достигая среза сопла они возмущают течение вниз по потоку, образуя таким образом замкнутый цикл с «накачиванием» интенсивности звукового давления в окрестности среза сопла. Отмечается, что этот процесс может привести к повреждению материала сопл. Прямое численное моделирование на основе уравнений Навье — Стокса и мультиоператорной схемы 7-го и 9-го порядков позволило описать этот эффект.
В расчетах использовалась более грубая сетка, чем обычно применяется в прямом численном моделировании (DNS), и мультиоператорные схемы 7-го порядка точности для получения поля течения плоской сверхзвуковой струи. Задача состояла в изучении акустических полей недорасширенной струи, в том числе и скрич-эффекта, при использовании обсуждаемого численного метода. Постановка задачи аналогична [20], где для ее решения использовался LES.
В расчетной области сопло моделировалось двумя параллельными пластинами. На входной границе задавались характеристические граничные условия вне сопла и фиксированные входные параметры потока с M = 1 и тонкими пограничными слоями. Заданное отношение давлений Pjet /= 2 соответствует числу Маха при полном расширении M = 1.52. Расчеты проводились при числах Рейнольдса 103 — 105. Re = 105 соответствует условиям эксперимента [21]. Представленные результаты получены на сетке с числом узлов 97x141. В результате расчета на начальном участке струи образуется система скачков («бочки»), которая с некоторого момента времени становится нестационарной (изолинии поля скорости показаны на рис. 13). При этом сдвиговые слои на границах струи после третьей от среза сопла «бочки» трансформируются в почти изолированные вихревые структуры. Подробный анализ этих структур, а также полей давления и плотности показал, что нестационарные ударные волны в третьей бочке «проникают» в зазоры
Рис. 13. Уровни поля скорости
Рис. 14. Поля модулей градиента плотности в моменты времени Г. а — t = 200; б — t = 202; в — t = 204; г — t = 206
между сохранившимися непрерывными сдвиговыми слоями и соседними вихревыми структурами, образуя волны сжатия вне струи. Эти волны, распространяясь вверх по потоку в виде полуокружностей, достигают среза сопла, влияя на течение за срезом. Они и являются скрич-волнами. На рис. 14 приведены поля модулей градиента в последовательные моменты времени (а — ґ = 200; б — ґ = 202; в — ґ = 204; г — ґ = 206). Они аналогичны шлирен-фотографиям, используемым в экспериментах. На рисунке можно видеть изменения с течением времени расположения сдвиговых слоев и ударных волн, а также акустических волн, вызывающих скрич-эффект.
Частота этих волн соответствует первому максимуму в спектре давления, показанному на рис. 15, а. Следующие максимумы соответствуют субгармоникам. Сравнение этого спектра со спектром, полученным с помощью ЬББ [20] при приблизительно аналогичных условиях показывает, что он имеет большее число локальных максимумов. Однако значение числа Струхаля, соответствующее частоте скрича, близко к значению, приводимому в [21].
На двух других частях рис. 15 частота скрича (б) и уровень звукового давления (в) сравниваются с экспериментом [21], теорией [22] и численным расчетом [20]. Видно, что полученная
б)
SPL, дБ
\ 170
\ 160 •
150 140 + • x^xxjt * + + + ** X * X + X * х +11
130
1.2
.4 1.6 М; 1.8
1
1.2
1/1
Рис. 15:
а — звуковое давление в зависимости от числа 8Ь; б — зависимость числа скрича от числа М; в — зависимость величины звукового давления от числа М;
□ □ □ — настоящий расчет; • • • — расчет [20]; х х х — эксперимент [21];-----------теория [22]
зависимость числа Струхаля от числа Маха хорошо согласуется с экспериментом. Сами же уровни звукового давления меньше, чем в эксперименте, согласно которому утолщение стенок сопла на срезе может привести к увеличению звукового давления на 20 дБ; в наших же расчетах стенки предполагались бесконечно тонкими.
ЗАКЛЮЧЕНИЕ
В работе представлена общая методика построения схем желаемого порядка с использованием линейных комбинаций базисных операторов (мультиоператоров). В силу образующихся при этом свободных параметров возникает возможность оптимизации схем с целью увеличения диапазона волновых чисел гармоник, разрешаемых сетками, для которого фазовые и амплитудные ошибки малы. Это делает подобные схемы привлекательными для задач вычислительной аэроакустики. Представленные результаты решения трех различных задач и их сравнения с экспериментальными данными иллюстрируют успешное применение мультиоператорных схем.
Оценивая их применимость в других задачах аэрогидродинамики, отметим, что схемы рассчитаны на использование равномерных сеток. Такие сетки могут учитывать сгущение узлов в физических областях только в случае достаточно гладких (например аналитических) преобразований координат. Поскольку гладкость решений является основным условием достижения высоких порядков любой методики, применение мультиоператорных схем в случае сквозного счета сильных разрывов становится проблематичным. Хотя они являются консервативными, возникающая немонотонность в окрестности скачков требует введения ограничителей потоков, понижающих точность решений. Вместе с тем, результаты решения задачи о скриче показывают, что при небольших числах Маха внутренний диссипативный механизм схем способен обеспечить монотонность решений.
Работа выполнена при поддержке РФФИ (грант № 08-01-00354-а) и проекта Программы фундаментальных исследований ОМН РАН № 3.
ЛИТЕРАТУРА
1. Shu C.-W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws // NASA CR-97-206253, ICASE Rep. N 97-65, Langley Res. Center, Hampton (1997).
2. Tam C. K. W., Webb J. C. Dispersion-relation preserving finite difference schemes for computational acoustics // J. of Comput. Physics. 199З. V. 107, p. 262 — 281.
3. Толстых А. И. Об итерационных схемах с нецентрированными компактными аппроксимациями // Докл. АН СССР. 1992. Т. З26. № З, с. 425 — 4З1.
4. Толстых А. И. Мультиопеpатоpные схемы пpоизвольного прядка, использующие нецентpиpованные компактные ап^оксимац^ // Докл. РАН. 1999. Т. З66, № З, с. З19 — З22.
5. Толстых А. И. О построении схем заданного порядка с линейными комбинациями операторов // Ж. вычисл. матем. и матем. физ. 2000. Т. 40, № 8, с. 1206 — 1220.
6. Липавский М. В., Т олстых А. И., Чигерев Е. Н. О схеме с мультиопера-торными аппроксимациями 9-го порядка для параллельных вычислений и о ее применении в задачах прямого численного моделирования // Ж. вычисл. матем. и матем. физ. 2006, Т. 46, № 8, с. 14ЗЗ — 1452.
7. Толстых А. И. О методе численного решения уpавнений Навье — Стокса сжимаемого газа в штоком диапазоне чисел Pейнольдса // Докл. АН СССР. 197З. Т. 210, № 2, c. 48 — 51.
8. Толстых А. И. Компактные pазностные схемы и их ^именете в задачах аэpо-гидpодинамики. — М.. Наука, 1990.
9. Толстых А. И. Об одном классе нецентpиpованных компактных pазностных схем пятого поpядка, основанных па аппpоксимациях Падэ // Докл. АН СССР. 1991. Т. З19, № 1, c. 72 — 77.
10. T o l s t y k h A. I. Development of arbitrary-order multioperators-based schemes for parallel calculations. 1. Higher-than-fifth-order approximations to convection terms // J. of Comput. Physics. 2007. V. 225, p. 2ЗЗЗ — 2З5З.
11. Толстых А. И. Об одном семействе компактных аппроксимаций и основанных на них мультиоператорных аппроксимациях заданного порядка // Докл. РАН. 2005. Т. 40З, № 2, с. 172 — 177.
12. T o l s t y k h A. I. Development of arbitrary-order multioperators-based schemes for parallel calculations. 2. Families of compact approximations with two-diagonal inversions and related multioperators // J. of Comput. Physics. 2008. V. 227, p. 2922 — 2940.
13. Tam C. K. W. Problem l-aliasing. fourth computational aeroacoustic (CAA) in. proceedings of the workshop on benchmark problems // NASA CP-2004-2159, 2004.
14. Monkewitz P. A., Bechert D. W., Barsikow B., Lehmann B. Self-excited oscillations and mixing in a heated round jet // J. of Fluid Mech. 1990. V. 21З, p. 611 — 6З9.
15. T am C. K. W., Burton D. E. Sound generated by instability waves of supersonic flows. Part 2 // J. of Fluid Mech. 1984. V. 1З8, Cambridge University Press, p. 27З — 295.
16. Kopiev V. F., Chernyshev S. A., Zaitsev M. Yu., Kuznetsov V. M. Experimental validation of instability wave theory for round supersonic jet, in Proceedings of The 12th AIAA CEAS Aeroacoustics Conference (27th AIAA Aeroacoustics Conference) 8 — 10 May 2006. — Cambridge, Massachusetts, American Institute of Aeronautics and Astronautics. 2006.
17. Troutt T. R.Mclaughlin D. K. Experiments on the flow and acoustic properties of a moderate-Reynolds-number supersonic jet // J. of Fluid Mech. 1982. V. 116, Cambridge University Press, p. 12З — 156.
18. Савельев А. Д. Расчеты течений вязкого газа па основе (д-у)-модели турбулентности // Ж. вычисл. матем. и матем. физ. 200З. Т. 4З. № 4, с. 589 — 600.
19. Lesshafft L., Huerre P., Sagaut P. Frequency selection in globally unstable round jets // Physics of Fluids. 2007. V. 19, 054108, American Institute of Physics, p. 1 — 10.
20.Berland J., Bogey C.Bailly C. Numerical study of screech generation in a planar supersonic jet // Physics of Fluids. 2007. V. 19, 075105, American Institute of Physics, p. 1 — 14.
21. Panda J., Raman G., Zaman K. B. M. Q. Underexpanded screeching jets from circular, rectangular, and elliptic nozzles // NASA TM-2004-212481, AIAA Paper 97-162З, 1997.
22. Tam C. K. W. The shock-cell structures and screech tone frequencies of rectangular and nonaxisymmetric supersonic jets // J. of Sound and Vibration. 1988. V. 121, Academic Press, p. 1З5 — 147.
Рукопись поступила 7/II2011 г.
