Обобщение метода Годунова на задачи вычислительной аэроакустики Текст научной статьи по специальности «Физика»

Том ХЫ

УЧЕНЫЕ ЗАПИСКИ ЦАГИ 2010

№ 1

УДК 534.83

ОБОБЩЕНИЕ МЕТОДА ГОДУНОВА НА ЗАДАЧИ ВЫЧИСЛИТЕЛЬНОЙ АЭРОАКУСТИКИ

И. С. МЕНЬШОВ

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

Ключевые слова: аэроакустика, численное моделирование, вариационная задача Римана.

В конце пятидесятых годов С. К. Годунов предложил оригинальный численный метод для задач газовой динамики, в основе которого лежит точное решение задачи Римана, более известной в механике как задача о распаде произвольного разрыва в газе [1]. С математической точки зрения она представляет собой задачу Коши для законов сохранения, определяющих движение сжимаемой жидкости, с кусочно-постоянным распределением начальных данных. С точки зрения физики задача Римана описывает взаимодействие двух сжимаемых потоков газа. Замечательной особенностью является тот факт, что она, в силу своей автомодельности, допускает точное аналитическое решение, представляемое в компактной явной форме.

Численная схема С. К. Годунова строится стандартным методом конечного объема. Возникающие при этом численные потоки на гранях счетных ячеек определяются значением дифференциального потока на точном решении соответствующей задачи Римана. Этот подход легко обобщается на многомерный случай [2]. Идеи С. К. Годунова получили очень широкое распространение во всем мире, и в настоящее время в литературе по вычислительной гидродинамике используется устоявшийся термин — методы годуновского типа [3].

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

1. Формулировка возмущенной задачи Римана. В классической постановке задача Римана сводится к задаче Коши для системы уравнений одномерной газодинамики:

^2+—=о, (1)

ді дх

где о — вектор консервативных переменных, Е — вектор потока с кусочно-постоянными начальными данными О = О/ при х < 0 и О = Ог при х > 0.

Эта задача имеет единственное решение при любых, но физически допустимых начальных данных, которое представляется кусочно-гладкой функцией автомодельного параметра X = х/ї и векторов О/ и Ог:

О(і, х) = (X, О/, Ог). (2)

Оно было получено Н. Е. Кочиным [4]. Детальное описание функции О можно найти в работе [2].

Возмущенная задача Римана формулируется следующим образом. Пусть известно решение

задачи Римана Ой (X, О/, Ог) при некоторых начальных данных О/ и Ог. Требуется найти, как

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

О/ ^ О/ +50/, Ог ^ Ог + 50г. (3)

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

= М/5О/ + М г 5ОГ, (4)

где матрицы М/ и М г, которые мы будем в дальнейшем называть вариационными, являются по сути якобианами функции решения задачи Римана по векторам начальных данных

М//г = М//г (X, )=^~. (5)

д°//г

Таким образом, возмущенная задача Римана сводится к нахождению вариационных матриц. Эти матрицы являются кусочно-гладкими функциями X, причем области гладкости у них те же, что и в основной задаче Римана, и определяются волновой структурой последней. Они не зависят

от X в тех подобластях, где в основной задаче реализуется однородный поток, т. е. в невозмущен-

ных областях и контактных зонах.

В зоне центрированной волны разрежения матрицы М^г определяются решением системы

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

1/г (X) = цвр ((X)), (6)

м

где цвр — вариационная матрица волны разрежения. Она приведена в Приложении 1.

Решение в контактной зоне определяется на основе соответствующих линеаризованных соотношений на ограничивающих эту зону разрывах. Тип разрыва зависит от конкретных начальных данных. Это может быть или слабый разрыв, или ударная волна. На слабом разрыве непрерывны энтропия и соответствующий римановский инвариант (точнее, его линейный аналог). На ударной волне уравнения для вариаций вытекают из линеаризации соотношений Рэнкина — Гю-гонио:

где А — якобиан функции потока К (Л = ЭК/ЭО), а квадратными скобками обозначается ска-

чок параметра на разрыве.

Анализ всех возможных волновых конфигураций, возникающих при распаде разрыва, показывает, что во всех случаях решение контактной зоны слева и справа от разрыва определяется с точностью до одной произвольной константы С:

Таким образом, чтобы завершить построение решения, необходимо определить константы Сі и Сг. Для этого имеются как раз два соотношения, выражающие непрерывность вариации давления и скорости на контактном разрыве. Окончательный результат можно представить в достаточно компактном виде. За деталями мы отсылаем читателя к работе [5].

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

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

Вводя обозначения «верхняя черта» и «крышка» для базовых величин и малых возмущений соответственно, вектор консервативных переменных представим в виде ц ц ц . Интегрируя уравнения Эйлера методом конечных объемов на некоторой заданной сетке и подставляя в полученные дискретные уравнения указанное выше разложение, можно прийти к следующей линейной модели, которая описывает распространение малых возмущений по заданному, вообще говоря, неоднородному, базовому полю:

где а — индекс грани счетной ячейки, юг- — ее объем, 5а — площадь грани, Та — матрица перехода к локальному базису. Стоящая в правой части этих уравнений переменная 8 включает все нелинейные члены. Они отвечают, главным образом, за генерацию малых возмущений и обычно моделируются специальным образом на основе экспериментальных данных или других расчетов. Этот вопрос требует специального рассмотрения и выходит за пределы настоящей работы; обсуждение его можно найти в [7].

Ключевым моментом в рассматриваемом методе является аппроксимация акустического

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

6ОІ/г = П1/г 8О//г + С1/г ,

(8)

где матрица и вектор т1г выписываются в явном виде для каждого типа разрыва (см. Приложение 2).

(9)

где Л (О ) = ЭЕ/ ЭО — якобиан невязкого потока, а О? = О? (, О?, О )) — решение базовой

задачи Римана на ребре с векторами О” и О?^), представляющими базовое течение на ребре со стороны текущей і и соседней ?(і) ячейки соответственно.

Л к

Вектор О? есть не что иное, как решение возмущенной задачи Римана, и записывается через вариационные матрицы:

О? = Мі (о, О?, і) )Г + МО(0 (о, о?, о°(і) )))

(11)

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

Интегрирование по времени получающейся в результате системы дифференциальных уравнений (9) выполняется по 3-уровневой схеме Рунге — Кутта с коэффициентами а0 = 1/3, а10 = 0, а11 = 2/3, в1 = 1/4, в2 = 0, в3 = 3/4, требующей минимального хранения данных и обеспечивающей при этом второй порядок аппроксимации [6].

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

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

Основное течение представляется двумя однородными потоками, связанными на фронте соотношениями Рэнкина — Гюгонио. Такое течение определяется числом М перед волной. Звуковое поле перед волной имеет вид:

где к — волновое число, ю = кс. Амплитуда падающей волны, нормализованная на значение давления перед скачком, составляла 0.01.

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

Результаты расчета (уравнения (9) с 8 = 0) представлены на рис. 1, а, где приведены численные значения акустического давления вдоль направления течения для волны с числом М = 1.5. Амплитуда звуковой волны увеличивается при переходе через скачок. При этом метод не дает каких-либо искусственных всплесков или провалов в распределении около фронта волны. Профили падающей и проходящей волн хорошо выдерживаются. Их амплитуда практически не меняется по мере продвижения вниз по потоку, хотя сетка была не столь мелкой и составляла 15 ячеек на волну.

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

0.02-

я 0.01

І -0.01 -

о

е

-0.02-

скачок ;

0-

-----------------і-------■--------1--------.-------І--------•--------г

0 2 4 6 8

а)

х

---1—і—1—і—1—і—<—і—1—і—1—і М

0 2 4 6 8 10 12

б)

Рис. 1. Прохождение звуковой волны через фронт ударной волны: а — распределение давления; б — коэффициент усиления в зависимости от числа М

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

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

источниковый член представляется суммой разностей потоков Е ( + ц) и Е +17 через все грани

счетной ячейки. На верхней и нижней границах ставятся периодические условия, на правой — граничные условия неотражающего типа [9]. Левая граница служит для генерации турбулентных возмущений. Сетка состоит из 512 X 256 равномерно распределенных ячеек.

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

В виртуальной области разыгрывается турбулентное поле возмущений, удовлетворяющее следующим свойствам: а) изотропность, б) Шу('V) = 0, в) р = р = 0, г) спектр кинетической

Аналогичная модель турбулентного поля используется во многих работах, связанных с моделированием турбулентных граничных условий. В работе [10] можно найти более детальное обсуждение этого вопроса.

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

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

энергии

(12)

терной скоростью турбулентных пульсаций у (или турбулентным числом Маха Мт = с).

Волновое число к определяет число волн, укладывающихся на линейном размере области Ь. Соответственно величина Ь/к определяет длину волны, а N/к — число ячеек на волну, где N — число ячеек расчетной области.

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

1) турбулентная кинетическая энергия е' = 0.5^ (V)

2) продольный тейлоровский масштаб Х'х = 2е^^^ ^ к0 ;

3) плотность вероятности турбулентной скорости Ро^1(V) и Рй¥2 (V), соответственно

в двух сечениях, расположенных перед фронтом волны (х = хск — 0.25Ь) и за фронтом (х = хск + 0.25Ь).

На рис. 2 показано распределение энергии пульсаций е' (а) и тейлоровского масштаба Х'х (б) в продольном направлении для значений параметров М = 5 и Мт = 0.08. Мелкие вихри дис-сипируют быстрее, чем крупные. Поэтому можно ожидать, что мелкомасштабная турбулентность (= = 60) будет затухать интенсивнее, чем турбулентность более крупных масштабов (к0 = 40 и

к0 = 20). Именно это и наблюдается в численных экспериментах. Энергия турбулентности с к0) = 60 падает почти на два порядка в сверхзвуковом потоке перед волной, в то время как для турбулентности к0 = 20 это падение составляет один порядок.

Ударная волна усиливает турбулентность. Турбулентная энергия скачком увеличивается. При этом крупные вихри оказываются более инертными по отношению к ударной волне, чем более мелкие. В силу этого обстоятельства скачок энергии на волне оказывается меньше в случае к0 = 20 по сравнению с к0 = 40 и к0 = 60, и интенсивность турбулентности во всех трех случаях за волной почти выравнивается (рис. 2, а).

Линейный масштаб турбулентности возрастает на ударной волне (рис. 2, б). Этот эффект более выражен в мелкомасштабной турбулентности, где он оказывается почти четырехкратным.

Характер турбулентности также меняется при переходе через ударную волну. На рис. 3 приводится вид плотности распределения вероятности для скорости в сечениях до и после волны. Видно, что в набегающем потоке турбулентность изотропна, тогда как за фронтом волны появля-

Рис. 2. Взаимодействие турбулентности с ударной волной: а — распределение турбулентной энергии; б — распределение линейного тейлоровского масштаба

Рис. 3. Взаимодействие турбулентности с ударной волной:

плотность распределения вероятности в набегающем потоке (слева) и за фронтом волны (справа) в случае крупномасштабной (а) и мелкомасштабной (б) турбулентности

Рис. 4. Мгновенное поле завихренности

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

Заключение. В статье изложен численный метод аэроакустики, который является прямым аналогом метода С. К. Годунова в газодинамике. В его основе лежит возмущенная задача Римана

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

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

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта 09-01-00851-а, 09-01-92102-яф-а).

Приложение 1. Вариационная матрица в области центрированной волны разрежения находится из линейной системы дифференциальных уравнений, являющейся следствием вариации соответствующих соотношений:

и ± а -Х = 0, си ± р = 0, / = 0,

(13)

где и, р, а, s — скорость, давление, скорость звука, энтропия, с = ра, штрих обозначает производную по автомодельной переменной Х = хД. Здесь и далее верхний знак следует брать для правой волны, а нижний — для левой.

Граничным условием для получающейся в результате системы будет служить условие сращивания на передней характеристике Х = Х0:

би + и'5Х0 = и0,

где вектор и = (и, р, 5). Решение этой задачи записывается в виде

би = Цврби0,

где матрица цвр имеет следующий вид [5]:

1 — Г ±(1 — Г)/с0 ±[ю(1 — Г) + Г(Эа/ Эх )р

Цвр =

±сГ

0

сГ/ с

0

сГ(- (да/ дх )р

(14)

и параметры Г и ю определяются формулами:

2

Г = -

7 + 1'

7=1+р

V

др

Ю= j с

Р0

-2 ( Эс Эх

(15)

В случае идеального совершенного газа параметр 7 совпадает с показателем адиабаты, а ю дается простым выражением:

ю=

а - а.

0

(7-1)

(16)

Приложение 2. Со ссылкой на работу [5] приведем формулы для матрицы N и вектора т, которые определяют решение в контактной зоне с точностью до одной произвольной константы: би = 0 + Ст. Выбирая в качестве этой константы вариацию скорости на контактном разрыве

би, из условия сращивания на слабом разрыве (задней характеристике волны разрежения) с [би ]±[бр ] = 0 и [бх] = 0 получим:

(17)

" 0 0 0 " " 1 "

N= ±с с/с0 сю , т= ±с

0 0 1 0

В случае ударной волны в качестве константы удобно выбрать комбинацию

где

где A тура;

C=

ôm0 +-^ m0, Ро

?0 = р0 (0 — Х0) = р( — Х0) — поток массы через фронт ударной волны. Тогда из условий сращивания (7) получим:

1 I1 + mQ )/m0 +X0 -m0a0T0 - X3 "-1-X1 "

N= 0 -m0 Xo m0 Хэ , m= X1m0

0 -A/(m0T ) T,/ T A/ T

(19)

i = u — М0 = mo/Co — относительное число Маха перед фронтом волны; T — темпера-

© = (3^ дp )s/T; параметры %15 %2> Х3 определяются выражениями:

хґ

О + m0aA

мО -1 :

X о =

1 + Mq + m0aA

0

m,

(мО -1)

Хэ m0T0

а + а

0

МО -1

(00)

ЛИТЕРАТУРА

1. Годунов С. К. Разностный метод численного расчета разрывных решений уравнений газовой динамики. — Мат. сборник. 1959. Т. 47. № 3, с. 271 —306.

2. Г одунов С. К., Забродин А. В, Иванов М. В., Крайко А. Н., Проко-повГ. П. Численное решение многомерных задач газовой динамики. — М.: Наука, 1976.

3. Godunov methods. Theory and applications / Под ред. E. F. Toro. — Kluwer Academic Plenum Publishers. — Manchester, UK. 2001, p. 1077.

4. КочинН. Е. Избранные труды. Т. 2. — М.: Изд. АН СССР, 1949, с. 5—42.

5. Men’shov I., Nakamura Y. On implicit Godunov's method with exactly linearized numerical flux // Computers & Fluids. 2000. Т. 29, N 6, p. 595—616.

6. Shuu C. W. Total-variation-diminishing time discretizations // SIAM J. Sci. Stat. Com-put. 1988. Т. 9, p. 1073 — 1084.

7. G o l d s t e i n M. E. On identifying the true sources of aerodynamic sound // J. Fluid Mech. 2005. V. 526. p. 337—347.

8. КонторовичВ. М. Отражение и преломление звука ударной волной // ЖЭТФ. 1957. Т. 33, с. 1527—1528.

9. Men’shov I., Nakamura Y. Implementation of the variational Riemann problem solution for calculating propagation of sound waves in nonuniform flow fields // J. of Comput. Phys. 2002. V. 182, p. 118—148.

10. Lee L., Lele S. K., and Moin P. Direct numerical simulation of isotropic turbulence interacting with a weak shock wave // J. of Fluid Mech. 1993.V. 251, p. 533—562.

Рукопись поступила 19/III2009 г.