Реализация алгоритма и результаты.

Введение

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

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

Численное решение задач аэроакустики основывается на уравнениях Навье-Стокса, Эйлера или их модификациях, продиктованных спецификой задачи. Целое семейство методов решений данных уравнений основано на усреднении, поэтому они становятся практически неприменимыми в некоторых задачах, так как, сглаживая, решения «стирают» информацию о высокочастотных волнах. Прямое численное моделирование остается не берущимся барьером для современных компьютеров, так как вычисления приходится выполнять на больших сетках с крайне мелким разрешением, что требует использования огромной памяти и скорости вычислений, пока технически недоступных. Собственно, поэтому, численное моделирование чаще применяется для калибровки более энергоэкономичных моделей, а не для решения прикладных инженерных задач.

Исходя из указанных выше трудностей, предлагается исследовать конкретно волновые аспекты задач. Внутренними источниками звука являются, очевидно, вихревые структуры и турбулентные возмущения среды. Однако, расчет звука, созданного турбулентностью, является чрезвычайно сложной задачей, не решенной до сих пор, несмотря на развивающиеся методы. Полное описание задачи распространения звуковых волны в турбулентной среде можно найти в [1, 2]. Мы же выпишем системы уравнений, описывающих данный процесс:

 

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

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

Вопрос моделирования полей стоит в науке достаточно давно, и первые практические работы начали появляться в 70-е годы [см. 2, 3], когда еще мощность компьютеров не позволяла без особых проблем получить практически полезные результаты. Число гармоник для расчетов в данных работах не превышало 100. В более поздних работах [4] их число увеличилось до 600, и дальнейшее увеличение числа гармоник не приводит к существенному увеличению точности. В данной работе были рассмотрены 10, 100 и 1000 бегущих волн.

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

Получения поля скоростей для свободной (вдали от твердых стенок) турбулентности является довольно сложной задачей, аналитическое или прямое численное решение которой довольно затруднительно. В общем случае данное поле является хаотическим, и даже явное задание такого поля невозможно. Вполне обоснованно можно рассматривать такое поля как стохастический набор бесконечного числа акустических волн разной частоты и амплитуды, но подчиняющихся некоторым статистическим условиям. Например, заданному спектру , где волновой вектор. В свободной турбулентности, ввиду отсутствия выделенных направлений, логично считать такой спектр изотропным, то есть , k = | k |. Также стоит отметить, что спектр во всех выражениях считается нормированным. Искомое поле пульсаций в трехмерном случае:

(1)

должно удовлетворять следующим условиям:

0) Бездивергентность. Это свойство является необходимым, если мы рассматриваем несжимаемую жидкость.

(2)

1) Поле должно быть однородным. То есть, среднее значение в каждой точке пространства и в любой момент времени должно быть равно нулю.

. (3)

Здесь и далее под средним значением понимается осреднение по реализациям. Выражение среднего значения скалярной или векторной величины , зависящей от параметров через интегрирование по реализациям, запишется как:

Где – соответствующие плотности распределений.

2) Должны быть выполнены соотношения:

а) Корреляционная матрица, в силу однородности и изотропичности поля, обязана иметь следующий вид:

где – функции из в . На диагонали корреляционной матрицы стоят функции одного и того же вида при фиксированных и . В частности, выполнено следующее:

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

Значения, стоящие

(6)

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

3) Должно быть выполнено спектральное соотношение:

(7)

где

есть прямое преобразование Фурье.

Описание алгоритма.

Искомое поле можно [1, 2] представить в виде разложения в ряд Фурье по гармоникам:

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

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

Возьмем дивергенцию от обеих частей этого тождества, и, с учетом свойства (2):

Откуда следует, что скалярное произведение , то есть прямое преобразование Фурье от искомого поля имеет вид

где – комплекснозначный вектор, который в общем случае зависит от волнового вектора и частоты. Его можно записать в виде суммы вещественной и мнимой частей:

(14)

Вектора и , естественно, имеют вещественные компоненты.

Поиск волновых векторов по заданному спектру не представляет особого труда, но нахождение амплитуд является более сложной задачей. Исходя из (2), (12) и (14) было предложено [1] искать амплитуды по формулам:

где под знаком «крест», как обычно, понимается векторное произведение, а определение векторов ξ и η было дано ранее. Нам хотелось бы, чтобы при N →∞ спектр нашего поля приближался к заданному. Благодаря такому заданию амплитуд, свойство бездивергентности поля (2) выполняется тождественно. В тензорных обозначениях соотношения (15) и (16) запишутся следующим образом (опустив индекс, обозначающий номер гармоники):

где под обозначен абсолютно антисимметричный тензор Леви-Чивиты.

Покажем, что поле, построенное с помощью данного алгоритма (9), удовлетворяет условиям (2)-(8). Для этого найдем первый и второй моменты поля (мат. ожидание и корреляционную матрицу).

Заметим, что среднее значение величины равно нулю.

Действительно, введем сферическую координатную систему, направляющим вектором которой является вектор . Тогда среднее значение, согласно (4), (10) и (17), запишется в следующем виде:

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

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

Выпишем некоторые свойства корреляционного тензора. В общем случае структурную функцию можно [1] для однородного изотропного поля представить в виде суммы компонент вдоль и поперек вектора r.

 

Для того чтобы найти дисперсию поля (корреляционную матрицу, посчитанную в нуле), для начала рассмотрим выражение вида (суммирование по повторяющимся индексам не производится, за исключением индексов в символе Леви-Чивиты):

/ m,n = 1,2,3 /

Теперь выделим отдельно произведения гармоник с одинаковыми номерами:

Вспомним, что символ Кронекера. Тогда первое слагаемое будет заведомо равно нулю в том случае, если . Теперь покажем, что среднее значение от второго слагаемого равно тоже нулю:

Усреднять выражение (20) нужно по всем случайным величинам, которые в это выражение входят: по и , причем, для каждого индекса и отдельно, так как разным индексам соответствуют разные случайные величины. Очевидно, среднее значение каждого из четыре слагаемых под знаком суммы в (22) равно в точности нулю, так как, усредняя по одному из чисел и , мы получим интегралы вида ( какая-то функция, не зависящая от ):

которые равны нулю каждый. Отсюда следует, что второе слагаемое в выражении (21) равно нулю, что и доказывает диагональную структуру ковариационной матрицы:

Теперь найдем значения , стоящие на главной диагонали этой матрицы. Для этого посчитаем:

Для краткости записи обозначим за фазу . Тогда (23) преобразуется в:

Очевидно, третье слагаемое в сумме будет равно нулю, и его можно не писать.

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

Теперь найдем среднее значение величины для каждого в сферических координатах:

Перекрестный член будет равен нулю.

(25)

где

Посчитаем сначала следующие интегралы тригонометрических функций:

Таким образом, интеграл (23) преобразуется в:

Нетрудно показать, что данные рассуждения приведут к таким же результатам для других компонент поля. И, таким образом, мы получаем, что среднее значение произведения двух одинаковых компонент поля:

Видим, что эта величина не зависит от радиус-вектора и времени в точке наблюдения.

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

(30)

Здесь мы сразу обнулили перекрестный член. Воспользуемся тригонометрическими тождествами:

Тогда выражение (28) запишется как:

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

Аналогично прошлым рассуждениям найдем величину (31) для первой компоненты. Очевидно, ее значение будет ровно таким же для второй и третьей компоненты. Запишем ее в интегральном виде, используя (4) и (10), и выбрав направляющий вектор сферической системы координат по вектору :

Посчитаем интегралы по компонентам и (перекрестный член равен нулю), затем по и по :

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

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

Реализация алгоритма и результаты.

На языке С++ была написана программа, со следующей функциональностью:

- получение случайных величин с заданным спектром.

- построение полей с заданным спектром.

- визуализация полученных полей.

- статистический анализ результатов.

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

Для тестов был выбран спектр:

Где за обозначен максимум спектра.

Число гармоник было выбиралось поочередно равным 10, 100 и 1000, расчетная сетка состояла из 40 делений по осям Х и по У, и 12 делений по Z. Визуализация представляет собой «фотографии» двумерных полей скорости на разных слоях по оси Z. Сетка равномерная, шаг сетки h был выбран равным 0.3, .

Моделирование случайных величин с заданной плотностью распределения выполнялось с помощью преобразования Смирнова (метода обратных преобразований). Моделирование гауссовых векторов производилось с помощью преобразования Бокса-Мюллера. Описание этих методов можно найти в [4]. Для машинной реализации равномерно распределенных случайных величин был выбран стандартный алгоритм получения псевдослучайных чисел двойной точности, вшитый в компилятор Visual Studio 2010.

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

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

Графики этих функций и их соотнесение с посчитанными в программе при различных числах гармоник можно найти на рис. 2(а, б, в) в Приложении.

Сравнение теоретического спектра с полученным в программе отражены на рис. 3.

Статистический анализ результатов показал, что, во-первых, при росте N (10, 100, 1000) спектр поля стремится к изначально заданному, а, во-вторых, корреляционные функции стремились к аналитически вычисленным.

Выводы.

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

2) Анализ полученных графиков показал, что при увеличении числа гармоник корреляционные функции стремились к аналитически вычисленным. Также, спектр стремился к теоретическому.

3) Было оценено примерное количество гармоник, при которых корреляционные функции становились едва различимы на графиках: N = 600~700.

Заключение.

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

Современные вычислительные возможности сполна позволяют получать практически применимые результаты. Одна реализация поля из 1000 гармоник занимает на персональном компьютере от 1 до 15 секунд, что позволяет проводить подобные вычислительные эксперименты, буквально, в домашних условиях.

Список использованной литературы.

[0] А. С. Монин, А. М. Яглом, Статистическая гидромеханика, М.: Наука. Физматгиз, ч.1,2,

[1] Ландау Л. Д., Лифшиц Е. М. Теоретическая физика. – Издание 5-е. – 2006. – Т. VI. Гидродинамика.

[2] Сабельфельд К.К., Решение краевых задач методом Монте-Карло. – Новосибирск: Наука, 1980, 154 с.

[3] Kraichnan R.H., Diffusion by a Random Velocity Field. – The Physics of Fluids, 1970, vol. 13, N 1, p. 22-31.

[4] R. Sandberg, H. Fasel, Direct Numerical Simulations of Transitional Supersonic Base Flows, AIAA Journal, Vol. 44, No. 4, April 2006.

[5] Боровская И.А., Моделирование случайных сигналов и полей в задачах аэроакустики. – Москва, Институт математического моделирования РАН, 2007, автореферат на соискание ученой степени кандидата физико-математических наук, 3-5 с.

[6] Вадзинский Р.Н. Справочник по вероятностным распределениям. - СПб.: Наука, 2001, 295 с.

[7] N. Metropolis, S. Ulam, The Monte Carlo Method, — J. Amer. statistical assoc. 1949 44 № 247 335—341.

[8] А. М. Обухов, Турбулентность и динамика атмосферы, Гидрометеоиздат, 412,[1] с. ил. 23 см Л. Гидрометеоиздат 1988.

Приложение.

Рис.1. Поля X и Y составляющих скоростей, построенные на различных высотах Z.

 

 

Рис.2. График теоретической поперечной и продольной составляющих корреляционной функции.

Рис. 2а. Продольная и две поперечных корреляция при числе гармоник = 10.

Рис. 2б. Корреляция при N = 100.

Рис. 2в. Корреляция при N = 1000 .

Рис. 3. Сравнение теоретического и посчитанного спектра при N = 100.

8. Контакты:

Стерлинг Григорий Григорьевич

Sterling239@mail.ru

8 (929) 503-52-08

 

Поделиться:

Воспользуйтесь поиском по сайту: