Моделирование акустических полей при обтекании тел потоком газа

Багапов Ариф Ренатович

Моделирование акустических полей при обтекании тел потоком газа

Бакалаврская работа

Математика

Дипломы

Вуз: Мордовский государственный университет имени Н.П. Огарева

ID: 6065e2c3e010f800016c945a

UUID: 89d4b450-752a-0139-ab50-0242ac18000a

Язык: Русский

Опубликовано: около 3 лет назад

Просмотры: 3

Багапов Ариф Ренатович

Источник: федеральное государственное бюджетное образовательное учреждение высшего образования "Национальный исследовательский Мордовский государственный университет им. Н. П. Огарёва"

Комментировать 0

Рецензировать 0

Скачать - 1,4 МБ

Поделиться работой

ФЕДЕРАЛЬНОЕ ЕОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ОБРАЗОВАНИЯ «НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ МОРДОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ИМ. Н. П. ОГАРЁВА» Факультет математики и информационных технологий Кафедра прикладной математики, дифференциальных уравнений и теоретической механики УТВЕРЖДАЮ Зав. кафедрой канд. физ.-мат. наук, доц. _______ Р. В. Жалнин 25 июня 2020 БАКАЛАВРСКАЯ РАБОТА МОДЕЛИРОВАНИЕ АКУСТИЧЕСКИХ ПОЛЕЙ ПРИ ОБТЕКАНИИ ТЕЛ ПОТОКОМ ГАЗА Автор бакалаврской работы 12.06.2020 А. Р. Багапов Обозначение бакалаврской работы БР-02069964-01.03.02-02-20 Направление 01.03.02 Прикладная математика и информатика Руководитель работы канд. физ.-мат. наук, доц. 12.06.2020 Р. В. Жалнин 18.06.2020 Д. К. Егорова Нормоконтролер канд. физ.-мат. наук, доц. Саранск 2020

ФЕДЕРАЛЬНОЕ ЕОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ОБРАЗОВАНИЯ «НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ МОРДОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ИМ. И. И. ОГАРЁВА» Факультет математики и информационных технологий Кафедра прикладной математики, дифференциальных уравнений и теоретической механики УТВЕРЖДАЮ Зав. кафедрой канд. физ.-мат. наук, доц. ________ Р. В. Жалнин 20 января 2020 ЗАДАНИЕ НА ВЫПУСКНУЮ КВАЛИФИКАЦИОННУЮ РАБОТУ Студент Багапов Ариф Ренатович, 403 группа 1 Тема Моделирование акустических полей при обтекании тел потоком газа Утверждена приказом № 286-с от 20.01.2020 2 Срок представления работы к защите 27.06.2020 3 Исходные данные для научного исследования (проектирования) Литература по газовой динамике, акустике, численным методам, программированию 4 Содержание выпускной квалификационной работы 4.1 Алгоритм моделирования течения газа 4.2 Расчёт акустических характеристик 4.3 Моделирование акустических полей при обтекании различных тел 5 Приложения 5.1 Приложение А (обязательное) Листинг программы Руководитель работы 20.01.2020 Р. В. Жалнин Задание принял к исполнению 20.01.2020 А. Р. Багапов 2

РЕФЕРАТ Бакалаврская работа содержит 62 страницы, 20 рисунков, 2 таблицу, 12 использованных источников, 1 приложение. ГАЗОВАЯ ДИНАМИКА, HLLC, ENO, WENO, PARAVIEW, ОБТЕКАНИЕ, АКУСТИЧЕСКОЕ ПОЛЕ, OASPL, КУБИЧЕСКИЙ ОБЪЕКТ, КАСКАД, СТУПЕНЬКА, КАВЕРНА. Цель работы - создать математическую модель акустического поля, возникающего при обтекании объектов и тел газом. При решении задачи использовалась среда разработки Microsoft Visual Studio 2019, язык программирования C++, пакет интерактивной визуализации Paraview 5.6.0, программа Microsoft Office Excel. Степень внедрения - частичная. Область применения - расчёты в области газовой динамики, акустики. Эффективность - данная работа направлена на изучение численных методов высокого порядка точности, а также на понимание процессов, происходящих при обтекании различных объектов потоками газа. 3

СОДЕРЖАНИЕ ВВЕДЕНИЕ 5 1. Алгоритм моделирования течения газа 6 1.1 Разностная схема для расчёта газодинамических течений 6 1.2 Вычисление дискретных потоков 8 1.3 Алгоритм реконструкции газодинамических параметров 12 1.3.1 ENO алгоритм 12 1.3.2 WENO алгоритм 16 1.4 Дискретизация по времени 20 2. Расчёт акустических характеристик 21 2.1 Анализ полученной модели с помощью пакета Paraview 21 2.2 Анализ колебаний давления 22 3. Моделирование акустических полей при обтеканииразличных тел 25 3.1 Один кубический блок 25 3.2 Каскад из двух кубических блоков 29 3.3 Ступенька 34 3.4 Каверна (выемка) 39 ЗАКЛЮЧЕНИЕ 44 СПИСОК ИСПОЛЬЗУЕМЫХ ИСТОЧНИКОВ 45 ПРИЛОЖЕНИЕ А (обязательное) Листинг программы 47 4

ВВЕДЕНИЕ В настоящее время методы математического моделирования повсеместно внедряются в производство, численный расчёт становится мощным способом повысить надёжность и точность инженерных исследований. Это касается, прежде всего, тех областей, в которых конструкторские ошибки стоят очень дорого: машиностроение, архитектура, освоение космоса и пр.. В частности, одной из сфер, в которых широко применяются численные методы моделирования, является авиация. Расчёт и анализ аэродинамических и акустических параметров воздушных судов полезен как в мирных, так и в военных целях. Расчёт данных параметров однозначно приводит к уравнениям газовой динамики, решение которых с помощью ЭВМ наиболее целесообразно в силу сложности и нелинейности этих уравнений. При решении практических задач зачастую приходится сталкиваться с газодинамическими течениями, характеризующимися нестационарностью, нелинейностью происходящих процессов, разнохарактерным и сложным механизмом взаимодействия, для моделирования которых необходимо использовать численные методы повышенного порядка точности, чтобы получить максимально близкие к реальным результатам. В данной работе рассматривается задача моделирования акустических полей при обтекании тел потоком газа. Для решения задачи использовались алгоритмы WENO для реконструкции газодинамических параметров, программный комплекс Paraview для сбора данных с полученной модели, программа MS Excel для анализа собранных данных при помощи формул, позволяющих интерпретировать колебания давления как акустические. 5

1. Алгоритм моделирования течения газа 1.1 Разностная схема для расчёта газодинамических течений Математическая модель изучаемых в работе процессов представляет собой закон сохранения масс, энергии и импульса [4]: др — + div(pv) = О, ^ + d iv n = 0 , (1.1.1) — + div(ve) = —div(pv) + Q + F ■v — div(lV), P = P(PJ ) , £ = £(P J ) ■ Здесь p - плотность газа, p - давление, - вектор скорости частицы, в - удельная внутренняя энергия на единицу массы, е = р I £ + — I - полная энергия, T - температура газа, П = pv0 v + рI - тензор плотности потока импульса (I - единичный тензор), Q - мощность тепловых источников, — * F - некоторая сила, действующая извне (например, сила тяжести), W - вектор плотности теплового потока. В качестве уравнения состояния принимается уравнение состояния идеального газа с показателем адиабаты у: 6

Р = (у- 1 )ер- ( 1 . 1 .2 ) В векторной форме система уравнений (1.1.1) без членов, которыми в исследовании пренебрегаем, приобретает вид: dU dt dF^U dF&U дх ду (1.1.3) где (1.1.4) При рассмотрении конкретной модели также необходимо задать начальные и граничные условия для полного описания решаемой задачи. Начальные условия показывают состояние системы в начальный момент времени, а граничные задают состояние системы на её границе: это могут быть условия непротекания, неприлипания, неограниченного втекания или вытекания. Для построения дискретной модели область непрерывного изменения аргумента накроем сеткой, равномерной по каждому направлению: Д — дх Х Ду, где (1.1.5) 7

Здесь Lx, Ly-размеры исследуемой области по соответствующим направлениям. Все газодинамические функции будем задавать как средние значения в ячейках сетки. Для аппроксимации системы (1.1.3) использовалась нелинейная консервативная квазимонотонная дифференциально-разностная схема [5]: H(l)л П н(2)-л— „(2) л /7 л — /7 /7 лиц I i+2j | Ц+2 V- = о , dt hУ ( 2) где , Нi7+1/2 (1.1.6) дискретные потоки, вычисляемые согласно пункту 1.2. 1.2 Вычисление дискретных потоков Дискретные потоки на границах между ячейками вычислялись с помощью схемы HLLC (Harten-Lax-van Leer-Contact). Впервые она была представлена Э. Торо в 1994 году. Данная схема является модификацией схемы HLL, заключающейся в учёте промежуточных волн. Основная идея обеих схем состоит в том, чтобы искать числовой поток напрямую, а не как в исходном методе Годунова, где сначала ищется решение Qt ±/ 2 для задачи Римана, которое используется для оценки интеграла по времени от потока и, отсюда, находится числовой поток. Для определения числового потока для двумерных уравнений Эйлера, который является вращательной нормальным инвариантности, для внутренности достаточно ячейки рассмотреть в силу расширенные одномерные уравнения Эйлера, выровненные в направлении нормали. Предположим, без ограничения общности, что нормальным является направление по оси Ox. Поэтому интересующие нас уравнения примут вид: 8

dQ dt где Q и F(Q) есть U и F Известно [6], dF(Q) _ (1.2.1) dx U из (1.1.4) соответственно. что вектор-столбец собственных значений соответствующего якобиана представляется следующим образом: ( и — с\ “ ) и + с/ (1.2.2) где с = I ■- скорость звука. В каждой ячейке мы рассматриваем локальную задачу Римана [11]: dcQ + axF{Q) = о, !?i = Q” , если х < 0, \ Qr = Q .,_, есл и х > 0. (1.2.3) Чтобы получить поток HLLC, рассмотрим волновую картину на рисунке 1, где присутствует промежуточная волна скорости S*. Рисунок 1 9

Предположим, что оценки скорости SL,S*,SR для трёх семейств волн нам известны (они будут описаны в работе ниже). Применение интегральной формы законов сохранения в области [лу,0]х[0,7] и [.xR,0]x[0,F] даёт, соответственно: F*l — Fl F*r — Fr + + Sl (Q*l ~ Ql )> S r (Q*r ~ Qr )> (1.2.4) где Q*K — \_P*K> P*KU>P*KV>P*Ke ]T> (1.2.5) где К = L, R. Это приводит к алгебраической системе с количеством неизвестных большим, чем количество уравнений. Один из способов решения этой проблемы состоит в том, чтобы ввести ряд предположений, согласующихся с точным решением задачи Римана (1.2.3), что приводит к следующим выражениям для векторов промежуточного состояния в (1.2.4): (1.2.6) где K = L, R. Тогда промежуточные потоки F*L и F*R в (1.2.4) полностью определены, и числовой поток имеет следующий вид: (1.2.7) 10

Для схемы HLLC необходимы оценки SL, SR и S*. Следующий выбор оценки, основанный на оценке давления в так называемой «звёздочной» области (области между SL и Sr) оказался полезным. Для начала необходимо найти оценки SL и SR: $L — UL ~ Sr — u r + cl 4 l > cr 4 r > (1.2.8) где 1, если p* < pK, _L Чк — (1.2.9) 1 + i f Gt _ l) ] 2 ' если p* > Pk - Здесь p* - оценка давления p* внутри «звёздочной» области, а K = L, R. Обратим внимание, что в (1.2.9) проводится различие между разрежением и ударом. То есть, если волна является волной разрежения (p* < pK), то скорость SK будет соответствовать характеристической скорости разрежения, которая является скоростью самой быстрого сигнала. Если имеется ударная волна , то скорость находится как приближение истинной скорости удара[10]. Что касается p*, то надёжным выбором является двухразреженное приближение: Это приближение получается, если априори предположить, что две волны давления (самые быстрые волны в решении задачи Римана) являются волнами разрежения. 11

После нахождения оценок SL и SR можно дать оценку S*: Pr ~P l +P lul (.sl ~ ul ) - P r ur (Sr ~ ur ) P l (s l ~ ul ) - P r (s r - (1.2.11) ur ) Как правило, предложенные оценки скорости были очень надёжными на протяжении многих лет [11]. 1.3 Алгоритм реконструкции газодинамических параметров 1.3.1 ENO алгоритм Алгоритм ENO неосциллирующие) (Essentially описывает Non-Oscillatory построение и - существенно локальный выбор аппроксимационной формулы с наименьшими осцилляционными свойствами (то есть наименьшей величиной колебаний) из набора формул того же порядка. Пусть задана равномерная сетка с шагом Xi <. 2 < X.I—1 <. 2 : (1.3.1.1) 2 2 будем считать, что эту сетку, при необходимости, мы можем продолжить вправо и влево, добавив элементы справа или слева соответственно. Пусть заданы средние значения некоторой функции v (х) в ячейках сетки: v ± <...< Vi <...< vN, где 12 (1.3.1.2)

Vi = Введём обозначение I ячейки ( xl_Ь/ 2 ,Xi+Ь/ 2) сетки. Для каждой ячейки I ( i = 1, iV) построим полином p(x) степени не большей К-1, который интерполирует функцию v (x) с порядком точности К в пределах данной ячейки, т.е. при р (x) = v (x) + О(AxK) , (1.3.1.3) ■Ь Д х JVx. + i} р (x) dx = v l~2 (1.3.1.4) и Тогда, обозначив v +_Ь/ 2 = р ( x l_ Ь/ 2) и v _+Ь/ 2 = р ( x l+Ь/ 2) , получим = v (x _ i ) + 0(AxK), l~ ) l_ \ v.+i = v {xi+i j + 0(AxK). (13.1.5) Такой полином получается [9], если рассмотреть первообразную V(x) функции v (x) : V(x) = ЛУ v (О df, и построить полином P(x) степени К, который будет интерполировать функцию V(x) на границах ячеек сетки. Положим р (x) = Р'( x) . 13

Представив полином P(x) в виде интерполяционного полинома Лагранжа [3]: к к X ^j-r+ f-l/2 а д = 2 ^ - г + т - 1 ) П ^ ~Г + 771—1 / 2 X;i-r+l-1/2 771= 0 1=1=0 0 1Фт при X = X;+х| 2 получим И г= оП q=0 (y i+ 1/2 £*0 P(^i+l/2) =Zf=01 | 2m=; +l x i - r + q - 1/2) q * m ,l Y\l= o{x i - r + m - l / 2 ~ x i - r + l - l / 2 ) Vi-r+j^Xi- Г+j ■ 1*0 Далее на равномерной сетке с шагом Дх для искомых значений заданной функции v (х) в точках х;_ ±| 2 и х ;+±| 2 получаем выражения: ^i—1/2 0 Crj^i-r+j > ^i+1/2 0 Crj^i-r+j > (1.3.1.6) где _ Y<K Т.*1=о U Kq=o ( r - q + 1 ) l* m q * m ,l Crj ~ L ™=j+1 n io m-l l*m Значения c;j для K = 0,1,2,...,5 приведены в таблице 1. 14 (1.3.1.7)

Т а б л и ц а 1 - Значения коэффициентов crj. K г J=0 j =1 J=2 J=3 J=4 -1/2 1/2 3/2 -7/6 5/6 5/6 -7/6 -23/12 13/12 7/12 -5/12 13/12 -163/60 77/60 9/20 -13/60 17/60 -21/20 1/3 -1/6 1/3 11/6 13/12 -5/12 7/12 13/12 -23/12 137/60 -43/60 47/60 47/60 -43/60 137/60 -1/4 1/12 -1/12 1/4 25/12 -21/20 17/60 -13/60 9/20 77/60 -163/60 1/5 -1/20 1/30 -1/20 1/5 137/60 1 -1 1 0 1 -1 0 1 -1 3 0 1 2 -1 0 4 1 2 3 -1 0 5 1 2 3 4 2 3/2 1/2 -1/2 11/6 1/3 -1/6 1/3 25/12 1/4 -1/12 1/12 -1/4 137/60 1/5 -1/20 1/30 -1/20 1/5 Таким способом можно получить K полиномов рг (х) ,г = 0,.. . ,К — 1, удовлетворяющих условиям (1.3.1.3), (1.3.1.4). Каждый полином будет соответствовать шаблону Sr Xi 1 2 ’Xi-r+K+jy (1.3.1.8) где г = 0,..., К —1 . Суть алгоритма ENO заключается в том, что из всех K возможных полиномов выбирается наиболее гладкий на отрезке [хt_±/ 2, х {+±/ 2] . Выбор осуществляется следующим образом [9]. Сначала рассматривается шаблон, состоящий из границ одной ячейки: 15

х i-1/2 X i+1/ 2 }» далее, если | ^ [xi-3/2 >x i- 1/ 2»xi+1/ 2 ] I < I^ [xi- 1/2» xi+1/ 2»x i+З/2 ] | > добавляется точка xi- 3/ 2, в противном случае xi+3/ 2. И так далее, пока в шаблон не будут включены K+1 точек. Здесь ^ [xi-l/2» ■■■»xi+s+l/2] — ^ [xi+l/2» ■■■»xi+s+l/2] ^ [xi-l/2» ■■■»xi+s-l/2] ) xi+s+l/2 —x i-1/2 ^ [xi- 1/ 2 ] —^ (xi- 1/ 2 ) . Такой способ выбора шаблона обусловлен представлением искомого полинома в виде интерполяционного полинома Ньютона: Р ( Х ) = K ( X j _ r _ 1 / 2 ) + K [X j_ r _ 1/ 2, X j _ r + 1 / 2 ] ( x - Xj_r _ 1 / 2 ) + ••• + V\x i-r -1/2» ■■■»xi+s+l/2] (x —x i-r - 1/ 2 ) ■■■(x —xi+s+l/2)Однако, более точные результаты даёт использование WENO алгоритма. 1.3.2 WENO алгоритм Алгоритм схемы WENO (Weighted ENO - взвешенные ENO) [9] построения интерполяционного полинома основан на ENO схеме, но в 16

отличие от неё использует для аппроксимации комбинации всех возможных шаблонов для данной ячейки: к- 1 (г) (1.3.2.1) (j)rv i+1/2’ v i+1/2 — I где к- 1 v il\/2 = 'Yj CiJVi-r+j> ;=о г = 0 ,, К 1. (1.3.2.2) Таким образом, необходимо найти весовые коэффициенты <нг, г = К — 1, удовлетворяющие условиям к- 1 (л)г > О> ^ г г=О 1, (1.3.2.3) и сохраняющие все свойства существенно неосциллирующих систем. Для гладкой на каждом шаблоне Sr, г = 0,.. .,К — 1, функции v (х) существуют константы такие, что к- 1 Vi+1/2 = 'Y j drVi+1/2 = v (xi+1/ 2 ) + 0(Лх2^-1). r=О (1.3.2.4) Например, для 1 < К < 3 имеем d0 = 1, k = 1; (1.3.2.5) d0 17

_ 3 _ 3 _ 1 _ dnо —~~7, d-1 — ~ , dn — к —3. 10' 1 5 ' ^ ~7~7, 10 Видно, что константы всегда положительны, а их сумма равна единице. Если (л)г = dr + 0(AxK -1),г = 0,..., К —1, (1.3.2.6) то получим, что кV £+1/2 = 1 1 (г) (x)rv i+1/2 _ г=0 v(x. i) + 0(Дх2*-1), г —0 , , К —1. (1.3.2.7) *+2 Действительно, ^ -i ^ -i ^ “>rVil\/2 - ^ г=О г=О tf-i = tf-1 = ^ <*>М1\/2 + v(xi+1/ 2) ^ ( 0)г - dr) - ^ drv{£(О + 1/2 г=О г=0 г=О = - ^г)(^+}1/2 - v(xi+1/2)) = ^ г=О _ 0(Ахк г)0{Ахк) = 0(Ах2К х). г=0 Если на каком-то шаблоне 5Го функция v (х) имеет разрыв, то для сохранения свойств ENO схем необходимо коэффициент а)Го сделать близким к 0. Поэтому весовые коэффициенты вычисляются следующим образом [8]: аг Шг —уК-1 < Zjs=о us r = 0 ,...,K - l, 18 (1.3.2.8)

где ar dy (s + (3r) 2 (1.3.2.9) Здесь (3r - так называемый «индикатор гладкости», s - малая положительная величина, введённая, чтобы избежать деления на 0. Далее полагаем s = 1 0 ""40 [7]. На (3r накладываются следующие условия: /Зг = О(Ах2), если функция v (х) является гладкой на шаблоне Sr, и Рт = 0 ( 1), если функция v (х) имеет разрыв на шаблоне Sr . Тогда получаем, что (л)г = 0(1), если функция гладкая и (л)г = О(Ах4), если функция имеет разрыв на шаблоне . Согласно [8] положим: Рг = д1рг (х) дх1 19 dx, г = 0,..., К —1, (1.3.2.10)

где pr (х) - полином, соответствующий шаблону Sr,r = 0,. ..,К — 1. Таким образом, для K = 3 имеем: 13.n9 1 _ n9 Ро = J 2 f a ” 2vi+i + v i+2 ) + 4 (3vi “ 4vi+ 1 + ^ + 2 ) f t = ^ ( f t- i - 2vf + v i+1) 2 + ^ (vf_i - v i+1) 2 , I3 _ f t = ^ f a “2 ~~ 2Vi~x + (1.3.2.11) 9 1 + 4 (3Vi~2 ~ f a - 1 + 3Vi) 1 При этом получаем WENO схему пятого порядка [8]. Положим для дальнейшего исследования K = 3. 1.4 Дискретизация по времени Для дискретизации по времени для уравнения вида: dU , N = ЮТ. (1-4.1) где L( U) - пространственный оператор, использовалась TVD-схема РунгеКутта третьего порядка [9]: U* = Un + At ■L{Un), U** = l u n +± U* +± A t - L ( U*), 1 2 Un + 1 = - U n +-U** + At- L(U**). 3 3 20 (1.4.2)

2. Расчёт акустических характеристик 2.1 Анализ полученной модели с помощью пакета Paraview Описанные выше схемы и подходы позволяют проводить моделирование течения газа около некоего объекта и, таким образом, получать плоскостно-временное распределение газодинамических переменных. Более всего нас будет интересовать давление и его пульсации, именно они и позволяют определить акустическое поле, как в ближней окрестности объекта, так и в любой точке моделируемой области. Наиболее важными характеристиками акустического поля с точки зрения разработчиков летательных аппаратов являются зависимость общего уровня звукового давления (Overall Sound Pressure Level - OASPL) от направления на точку наблюдения и спектральный состав акустического сигнала в точках наблюдения. Спектральные характеристики сигналов могут быть получены с помощью методов, основанных на дискретном преобразовании Фурье. В данной же работе мы подробно остановимся на расчёте общего уровня звукового давления. Введём дискретную сетку для времени: o)At = {A0 i = 0, - ti_lf |Aj| = At,AtNt = T], (2.1.1) здесь [0,7] - промежуток времени, на котором проводится исследование. Отсюда, пульсации давления р '( R, t) в точке R для каждого момента времени находятся по формуле: р'(Д, tj) = p(tf) - р(и _ г), i = 1,... Ntl 21 ( 2 . 1. 2)

Для получения колебаний давления р (tj) на отрезке времени [0,1] воспользуемся пакетом анализа и визуализации данных Paraview 5.6.0, и, конкретно, инструментом “Plot Selection Over Time”, который позволяет построить график изменения всех используемых переменных в исследуемой точке. Далее нам необходимо получить численные данные о давлении в исследуемой точке. Для этого, выбрав полученный график, перейдём к виду в форме таблицы (“Spreadsheet view”) и, отметив нужные столбцы (время и давление), экспортируем их в формат *.csv для удобства последующего расчёта в MS Excel. Аналогичную операцию повторим для некоторого количества точек, равномерно расположенных вокруг исследуемого объекта на равном расстоянии от его центра, чтобы в дальнейшем получить диаграмму направленности общего уровня звукового давления. Увеличение количества точек для рассмотрения сделает диаграмму более подробной, но вместе с этим увеличится и время расчёта необходимых параметров. 2.2 Анализ колебаний давления Полученный акустический сигнал (последовательность пульсаций давления) для каждой исследуемой точки подвергается дальнейшему анализу для получения интересующих акустических характеристик. Одной из наиболее показательных характеристик является общий уровень пульсаций давления, измеряемый в децибелах. Данный показатель определяется по следующей формуле [12]: ( 2 . 2 . 1) 22

где р о —2 -1 О 5 Па - минимальный порог слышимости звука, (р '( R) 2) - среднее значение квадрата колебаний давления в точке R, вычисляемое по формуле: <Р '(R ) 2) = С точки зрения инженерных (2.2.2) разработок наиболее полезным представлением о картине общего уровня звукового давления является азимутальная диаграмма направленности - график в полярных координатах, на котором отражена зависимость общего уровня звукового давления в некоторой плоскости для некоторого диапазона значений азимутального угла. Как правило, значения угла берутся либо от 0 до л, либо от 0 до 2л . В данной работе исследуется двумерный случай, поэтому выбор плоскости представления ограничен одной. Для данного исследования вычисление средних значений давления, общего уровня пульсаций давления, а также визуализация полученных данных осуществляется при помощи MS Excel. Строго говоря, в отличие от аэродинамических характеристик течения, которые довольно точно предсказываются и на малом масштабе, для расчёта шума необходимо использовать расчётную область, размер которой в поперечном направлении совпадает с размером экспериментальной модели. Однако в большинстве случаев это требует очень больших вычислительных ресурсов и в реальности оказывается практически неосуществимым. В связи с этим для расчёта шума, создаваемого при обтекании «квазидвумерных» тел на основе аэродинамических расчётов, выполненных в области исследования размером L, значительно меньшей, чем реальные размеры рассматриваемого в эксперименте тела, предложен ряд специальных поправок. 23

Простейшая из них базируется на предположении о том, что акустические сигналы от отдельных отрезков тела длиной l, на которые может быть разбито реальное тело, являются некоррелированными между собой. В этом случае суммарная спектральная мощность шума равна сумме спектральных мощностей от отдельных источников длиной l, то есть искомая поправка определяется соотношением (2.2.3) А(дБ) = 101g(iV), где А(дб) - величина поправки в децибелах, ехр где lexp - размер тела в эксперименте [1]. Например, для реального кубического тела со стороной lexp = 10 см, которое в расчётах представлено моделью со стороной l = 0.1 см искомая поправка будет равна Как можно понять, это достаточно грубое предположение, которое может привести к существенному занижению шума. Существуют и более точные поправки. Они, тем не менее, в данной работе не описываются из-за относительной трудности расчёта при использовании описанного выше метода анализа. 24

3. Моделирование акустических полей при обтекании различных тел 3.1 Один кубический блок Основой для расчётов стала программа, написанная на языке C++, листинг которой приводится в приложении А. В программе задаются начальные и граничные условия, параметры моделируемой области, размер и положение обтекаемых тел, а также условия на их границах. Результаты расчётов выводятся в формате *.vtk и позволяют визуализировать и наблюдать в программном комплексе Paraview все переменные, описанные законами сохранения на протяжении всего времени моделирования. Практическую часть исследования было решено начать с базового примера обтекания одного блока простой формы. Для моделирования данного процесса и всех следующих были выбраны следующие параметры: - размеры моделируемой области: 1,78D x 6,5D, где D - некоторый множитель для масштабирования, равный для данной модели 0,05 м. На данной области вводится сетка размерами 89х325. - газодинамические параметры: начальное состояние области: -5 р = 1, 20 кг/м - плотность; (u; v) = (0; 0) м/с - вектор скорости; p = 101325 Па - давление; параметры втекающего газа: -5 р = 1, 65 кг/м - плотность; (u; v) = (0; 114,4) м/с - вектор скорости; p = 158900 Па - давление; - шаг по времени: т = 2.0e-7с; 25

- размеры объекта в примере: 10х10 ячеек сетки, что соответствует 0,2D*0,2D (сторона равна 0,01 м). Так как сравнения с реальным проведённым экспериментом в работе не проводится, то невозможно внести корректную поправку в соответствии с формулой (2.2.3). Поэтому, считая, что размер экспериментального тела соответствует размеру моделируемого, поправку вносить не следует. По окончании расчётов модель была визуализирована в Paraview. Отдельные этапы протекания процесса показаны на рисунке 2. а) t = 0,0016 c в) t = 0,0058 c б) t = 0,0029 c Рисунок 2 В момент t « 0,0029 с, как видно на рисунке 2.б, происходит смещение области пониженного давления, которая до этого момента была симметричной (рисунок 2.а), вследствие превышения критического значения числа Рейнольдса, из-за чего в свою очередь образуется турбулентное течение. Это приводит к образованию так называемой вихревой дорожки (дорожки Кармана), которая видна на рисунке 2.в. Более явно её можно наблюдать на рисунке 3, где показано поле плотности исследуемой области. 26

Это явление для плохообтекаемых тел возникает при ограниченных значениях числа Рейнольдса и бывает довольно опасным для различных объектов (зданий, сооружений, техники), так как приводит к их разрушению. Поэтому важным моментом создания конструкций является учёт возможного появления вихревых дорожек и принятия мер по предотвращению этого. Рисунок 3 - Визуализация вихревой дорожки на графике плотности Далее, как было описано в главе 2, был использован инструмент “Plot Selection Over Time”, чтобы получить данные о колебаниях давления в точках, окружающих объект. Опытным путём было выяснено, что для данной модели оптимальным выбором являются 12 точек исследования, расположенных на равном расстоянии друг от друга. Также, для исследования влияния дальности расположения точки от объекта на уровень звуковых колебаний, было проведено несколько вариантов расчётов, в которых точки располагались на расстоянии 0,3Д 0,5D и 0,8D от центра объекта. 27

Рисунок 4 - Расположение точек исследования Полученные данные в численном формате экспортировались в MS Excel, и были проведены расчёты для определения общего уровня звуковых пульсаций давления. В итоге получены конечные данные, представленные в таблице 2, а также в виде азимутальной диаграммы направленности на рисунке 5. Т а б л и ц а 2 - Зависимость OASPL от направления на точку и расстояния до неё Азимут 180 150 120 90 60 30 0 330 300 270 240 210 Расстояние до центра объекта 0.3D 0.5D 0.8D 119,8097 127,9251 130,1272 129,8959 130,6925 131,4455 130,5572 130,0003 131,741 130,9154 130,6243 130,4161 131,6899 130,8593 130,2191 132,451 131,7413 131,2284 132,5678 132,2297 131,8027 131,7382 132,494 132,449 130,8002 132,5274 132,391 130,147 131,2535 131,2209 130,191 130,6873 130,7772 129,4627 129,5303 130,8049 28

Рисунок 5 - Азимутальная диаграмма направленности общего уровня звукового давления По приведённым таблице и диаграмме можно сделать вывод, что исследование на близком расстоянии более информативно, так как оно показывает уменьшение общего уровня шума за объектом, как следствие появления там области пониженного давления без последующих колебаний (до образования завихрений). Для близкого поля течения газа около объекта следует использовать более подробные математические модели, что и было проделано в данной работе, так как для этой области существенны нелинейные эффекты и могут образовываться скачки, где важна динамика вихреобразования, учет турбулентности и т.п. В свою очередь в дальнем поле, где подразумевается однородное течение вдали от источника шума, перенос акустических возмущений можно описывать волновым уравнением. 3.2 Каскад из двух кубических блоков Настоящий пункт практической части работы посвящён анализу тех же параметров, что и предыдущий, для каскада из двух объектов размером 0,2D x 0,2D, находящихся на расстоянии 0,3D друг от друга. На рисунке 5 приведены отдельные моменты визуализации процесса. 29

a) t = 0.0016 c б) t = 0.0029 c в) t = 0.0058 c Рисунок 6 На начальном этапе две области пониженного давления симметричны относительно продольной оси пары объектов, периодически объединяются в одну. Однако за вторым объектом область гораздо меньше в силу того, что он находится в области действия так называемого «кармана» за первым объектом, следовательно, испытывает меньшие нагрузки со стороны среды. Подобный эффект «кармана» часто используется автогонщиками на трассах и обычными автомобилистами, когда они пристраиваются за едущим впереди автомобилем. Это позволяет поддерживать скорость и экономить топливо за счёт того, что давление впереди автомобиля оказывается ниже, чем сзади. Спустя время, как и в случае с одним объектом, наблюдались вихревые дорожки Кармана (подробнее показаны на рисунке 7). В данном случае оба объекта будут под угрозой, так как вихри, срываясь с первого, наталкиваются на второй, который в свою очередь также порождает вихревые потоки. 30

Рисунок 7 - Визуализация вихревой дорожки на графике плотности Для расчёта акустического поля целесообразно использовать не окружность, а вытянутую фигуру, чтобы захватить оба объекта. Например, можно использовать эллипс, фокусами которого являются центры объектов. Количество точек исследования было увеличено до 16. Рисунок 8 - Расположение точек исследования 31

На рисунке 9 представлена полученная для каскада кубов диаграмма направленности общего уровня шума (OASPL). Рисунок 9 Как можно судить по диаграмме, минимальный уровень звуковых колебаний наблюдается сразу за вторым объектом, а также ближе к центру пространства между объектами, что обосновано наличием там областей низкого давления на протяжении всего времени. Заметим, что диаграмма не совсем симметрична: связано это с тем, что образование вихревых дорожек происходит преимущественно на правой стороне из-за чего в одних и тех же областях наблюдается как высокое, так и низкое давление, что в свою очередь делает картину менее информативной на больших промежутках времени, сглаживая средние значения. С левой же стороны ситуация более стабильна, а потому об уровне шума можно судить даже по расчётам приведённым выше. Отдельный интерес представляет сравнение графиков изменения давления за каждым из объектов, которое представлено на рисунке 10. По этим графикам можно видеть, как в самом начале временного отрезка ударная волна гасится первым кубом и оказывает на второй меньшее давление. Далее видно, как до момента (t = 0,0003 c) давление в обеих точках 32

постепенно падает с незначительными колебаниями. После того, как происходит смещение области пониженного давления и образуются вихревые потоки, колебания давления за вторым блоком становятся сильнее, что обусловлено срывами с этого блока сильных вихрей. Говоря об этих процессах в терминах акустики, ударная волна повлечёт за собой повышение уровня шума до пика, далее он будет спадать в обеих областях. После появления завихрений колебания снова начнут вызывать повышение уровня шума, но на данном отрезке времени между блоками это будет заметно меньше, так как там сохраняется область низкого давления. Тем не менее, уровень шума незначительно повысится и в этой области, что и показывает увеличение колебаний в конце временного отрезка. 0 0,001 0,002 0,003 0,004 0,005 0,006 Время, с Точка за первым объектом (м еж ду объектам и) Точка за вторым объектом Рисунок 10 - Графики колебания давления Расчёт общего уровня звуковых колебаний для точки между объектами дал значение 128,6 дБ - такая незначительная разница со значением для точки за вторым объектом (128,7 дБ), несмотря на меньшие колебания, обусловлен большим пиком в начале. 33

3.3 Ступенька Моделирование обтекания прямой ступеньки, как и расчёт создаваемого при этом акустического поля, является одной из важных задач аэродинамики. В том или ином виде эта базовая задача встречается во многих реальных проектах. Высота подъёма ступеньки, которая была выбран для модели, равен 0,5D. Подъём находится на расстоянии 2,4D от границы, через которую втекает газ. Обтекание подобных тел вызывает колебания давления и звуковые волны, как в ближнем, так и в дальнем поле, но в данной работе акцент сделан именно на моделирование ближнего поля. На рисунке 11 представлены этапы визуализации модели в программном комплексе Paraview. Можно наблюдать, как образуется отрыв потока и область пониженного давления, которая в дальнейшем увеличивается и сдвигается от края ступеньки. - 1,9е+05 - 180000 160000 1 140000 120000 - 100000 - 8.4е+04 а) t = 0.0003 c б) t = 0.002 c Рисунок 11 34 в) t = 0.006 c

Точки исследования были установлены, как показано на рисунке 12. Это расположение обусловлено желанием рассчитать общий уровень шума в непосредственной близости к ступеньке. Рисунок 12 - Расположение точек исследования На рисунке 13 представлен график зависимости общего уровня звукового давления от положения точки. Рисунок 13 - График зависимости общего уровня звукового давления от положения точки 35

Ударная волна, пришедшаяся на ступеньку, вызвала колебания давления, которое привело к большому уровню звукового давления. Низкий уровень шума в точке 8 объясняется тем, что возникшая на подъёме ступеньки область низкого давления постепенно увеличивалась и медленно двигалась в направлении потока. Также при контакте потока со ступенькой происходит отрыв, из-за которого основная часть потока проходит над поверхностью ступеньки, не касаясь её и не вызывая колебаний давления. Продемонстрировать разницу уровней звукового давления до ступеньки, на подъёме и после него позволит расчёт акустического поля на поверхности правой стенки, а также на подъёме ступеньки. Схема расположения точек показана на рисунке 14. Рисунок 14 - Расположение точек исследования В результате был получен график зависимости общего уровня звукового давления от расположения точки (рисунок 15), а также графики колебания давления в точках 1, 4, 6, 7 и 10 (Рисунок 16), чтобы подробнее разобрать происходящие в этих областях процессы. 36

Точки исследования Рисунок 15 - График зависимости общего уровня звукового давления от положения точки 250 230 210 £ 190 * 1 Ф 4 1 170 6 Ф 5 7 150 10 130 110 90 0 0,001 0,002 0,003 Время, с 0,004 0,005 0,006 Рисунок 16 - Графики колебания давления Видно, что общий уровень звукового давления согласуется с графиком, полученным для ближнего поля около ступеньки, то есть он снижается по направлению потока после подъёма. Довольно резкое возрастание уровня 37

шума в точке 10 показывает, что за моделируемый отрезок времени в данной точке область высокого давления сменилась зоной низкого давления, в то время как точки 7, 8, 9 находились в более стабильной части, где давление в основном понижалось. На рисунке 16 можно увидеть некоторые закономерности поведения давления в разных областях ступеньки. Области до подъёма ступеньки и после него характеризуются своего рода семействами графиков: можно увидеть, насколько похожи графики разных краёв одной зоны (1 и 4 и 7 и 10). Разница графиков заключается в следующем: - до подъёма: в общем уровне и пиковом значении давления, форма колебаний же практически одинакова. Чем ближе точка к подъёму, тем выше в ней уровень давления. Это и объясняет постепенный пологий подъём общего уровня звуковых колебаний на рисунке 15. - после подъёма: в размере так называемого «провала». Влияние на график в данной области оказывает расширяющаяся и двигающаяся область низкого давления. Чем дальше точка от края ступеньки, тем шире будет «провал». Точка 6 находится на выступе подъёма, именно там зарождается область низкого давления - это объясняет небольшой подъём и резкий «провал» графика в этой точке. Далее давление здесь возрастает и довольно сильно колеблется, что становится причиной высокого уровня акустического излучения. До точки 7 область низкого давления доходит спустя время и задерживается дольше, вызывая похожие эффекты. Точка 10 характеризуется начальным колебанием давления, но спустя время давление там лишь снижается, и этот процесс продолжается до завершения моделирования. 38

3.4 Каверна (выемка) Нестационарное пульсирующее течение, возникающее при обтекании плоской каверны параллельным ей набегающим потоком интересно как с точки зрения аэродинамических нагрузок, так и сильного звукового излучения. Несмотря на простую геометрию области, в ней имеет место довольно сложное течение, включающее турбулентные слои, возвратно циркуляционную зону, вихреобразование. Практическая ценность исследования этих процессов, заключается в том, что такую структуру, как каверна, можно использовать как средство управления сложными течениями, например в качестве стабилизатора пламени [2]. Размеры моделируемой каверны составляют 0,5D x 1D, её левый край находится на расстоянии 2,4D от границы, через которую втекает газ. В ряде изученных при подготовке работ рассматривается влияние размеров и формы каверны на создаваемое излучение, и оно, безоговорочно, существенно, однако было решено остановиться на простом примере прямоугольной каверны. Также существует ряд статей, в которых подробно описано, как пульсирующие потоки газа создают акустическое излучение, распространяющееся в дальнее поле, отражаясь от стенок каверны. В данной же работе будет рассмотрено колебание давления в ближнем поле каверны, а также внутри неё. На рисунке 17 представлен фрагмент визуализации модели. Рассмотрим его подробнее, чтобы понять, какие процессы здесь имеют место. На данном фрагменте поток движется слева направо. Внутри каверны можно наблюдать систему волн давления. Одна зона повышенного давления формируется у задней (правой) стенки, другая приближается к передней (левой). Между ними можно видеть две небольшие зоны пониженного давления и одну зону большего размера, которая опустилась в каверну вследствие разбиения начального вихря, образованного ударной волной, на 39

две: вторая часть вихря преодолела каверну. Сразу же за каверной наблюдается зона пониженного давления, вызванная расширением газа при преодолении им задней стенки. Чуть дальше по направлению потока можно заметить следы предыдущего распавшегося вихря [2]. - 1,6е+05 - 150000 Ф - 140000 | Ф а: - 130000 - 120000 - 1.1е+05 Рисунок 17 Как и в случае со ступенькой, для каверны будут проведены расчёты в ближнем поле над ней и на поверхности: до передней стенки, на ней, на дне, на задней стенке и после неё. Расположения точек исследования в обоих случаях показаны на рисунке 18. Зона пониженного давления, находящаяся в глубине каверны, остаётся там на протяжении всего времени моделирования, вращаясь по часовой стрелке, что иллюстрирует создаваемые в каверне циркуляционные потоки. 40

а) б) Рисунок 18 - Расположение точек исследования На полученном графике (рисунок 19) для точек над каверной видно, что общий уровень звукового давления в различных точках ближнего поля отличается не более чем на 0,7 дБ, что является пренебрежимо малой величиной по сравнению с абсолютными значениями. Следовательно, можно считать, что в ближнем поле наблюдается однородное постоянное звуковое излучение, распространяющееся в дальнее поле. Точки исследования Рисунок 19 - График зависимости общего уровня звукового давления от положения точки 41

Фактически же акустические волны, исходящие от каверны, условно делятся на 2 группы: первая распространяется от левой кромки и генерируется потоком давления внутри каверны, вторая - от задней (правой) кромки и возникает из-за взаимодействия вихревых потоков с этой кромкой. Разница между этими группами практически не заметна (на рисунке 17 можно обратить внимание на две пересекающиеся дуги волн, отходящих от кромок) на выбранном масштабе и области исследования, так как они сливаются в одну волну, уравниваясь при расчёте общего уровня звукового давления. Точки исследования Рисунок 20 - График зависимости общего уровня звукового давления от положения точки На рисунке 20, который показывает зависимость общего уровня звукового давления от положения точек на поверхности каверны, видно, как уровень звукового давления слегка поднимается, подходя к кромкам каверны. В промежутке между точками 2 и 10 можно наблюдать влияние на общий уровень звукового давления нелинейных эффектов, являющихся 42

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

ЗАКЛЮЧЕНИЕ В ходе выполнения выпускной квалификационной работы была разработана программа на языке C++, позволяющая моделировать обтекание различных тел потоком газа на основе законов сохранения. Использовались схемы реконструкции газодинамических параметров и расчета потоков высокой точности. Преимуществом созданной программы является возможность простого изменения параметров течения газа, параметров среды, размера и положения обтекаемых объектов. Был исследован способ моделирования акустического поля, создаваемого при обтекании объекта, при помощи расчёта общего уровня звукового давления(ОАБРЬ). Было проведено моделирование процесса обтекания и расчёт акустического поля для различных тел и объектов: куба, каскада из двух кубов, ступеньки, каверны. Для каждого объекта был дан анализ происходящих процессов на основании визуализации и графического представления акустического поля. Для куба было проведено исследование влияния удаления точки наблюдения акустических эффектов в пределах исследуемой области, которое показало, что изучение ближнего поля в приведённой области исследования наиболее информативно. 44

СПИСОК ИСПОЛЬЗУЕМЫХ ИСТОЧНИКОВ 1. Гарбарук А. В. Расчет аэродинамики и шума при обтекании тандема цилиндров / А. В. Гарбарук, Ф. Р. Спаларт, М. Х. Стрелец, М. Л. Шур // Математическое моделирование. - 2014. - Т. 26, № 6. С. 119-136. 2. Савельев А. Д. Численное моделирование акустического излучения двумерной каверны в дозвуковом потоке / А. Д. Савельев // Ученые записки ЦАГИ. - 2014. -Т. 45, № 1. - С. 57-74. 3. Самарский А. А. Численные методы : учеб. пособие для вузов / A. А. Самарский, А. В. Гулин. М. : Наука, 1989. - 342 с. 4. Самарский А. А. Разностные методы решения задач газовой динамики : учеб. пособие для вузов / А. А. Самарский, Ю. П. Попов. М. : Наука, 1980. - 424 с. 5. Тишкин В. Ф. Разностные схемы трехмерной газовой динамики для задачи о развитии неустойчивости Рихтмаера-Мешкова / В. Ф. Тишкин, B. В. Никишин, И. В. Попов, А. П. Фаворский // Математическое моделирование. - 1995. - Т. 7, № 5. - С. 15-25. 6. Courant R. On the solution of nonlinear hyperbolic differential equations by finite differences / R. Courant, E. Isaacson, M. Rees // Communications on Pure Applied Mathematics. - 1952. - Vol. 5, № 3. P. 243-255. 7. Henrik A. K. Mapped weighted essentially non-oscillatory schemes: Achieving optimal order near critical points / A. K. Henrik, T. D. Aslam, J. M. Powers // Journal of Computational Physics. - 2005. - Vol. 207, № 2. P. 542-567. 45

8. Jiang G. -S. Efficient implementation of weighted ENO schemes / G. -S. Jiang, C. -W. Shu // Journal of computational physics. - 1996. - Vol. 126, № 1. - P. 202-228. 9. Shu C. -W. Essentially non-oscillatory and weighted essentially non oscillatory schemes for hyperbolic conservation laws / C.-W. Shu // ICASE Report 97-65. - 1997. - 84 p. 10. Toro E. F. Riemann Solvers and Numerical Methods for Fluid Dynamics : A Practical Introduction / E. F.Toro. - 3rd edn. - Heidelberg : Springer-Verlag Berlin Heidelberg, 2009. - 724 p. 11. Toro E. F. The HLLC Riemann solver / E. F. Toro // Shock Waves. - 2019. - №29. - P. 1065-1082. 12. Wijker J. J. Spacecraft Structures / J. J. Wijker. - 1st edn. - Heidelberg : Springer-Verlag Berlin Heidelberg, 2008. - 504 p. 46

ПРИЛОЖЕНИЕ А (обязательное) Листинг программы #define _CRT_SECURE_NO_WARNINGS #include <cstdio> #include <cstdlib> #include <cstring> #include <math.h> double const M_PI = 3.14159265358979323846; const const const const const const const const const const const const const double D = 50.e-3; int STEP_OUT = 100; int NX = 89; int NY = 325; double XMIN = 0.0; double XMAX = 1.78 * D; double YMIN = 0.0; double YMAX = 6.5 * D; double HX = (XMAX - XMIN) / NX; double HY = (y MAX - YMIn )/ NY; double TMAX = 4.5e0; double TAU = 2.0e-7; int K_WENO = 3; const const const const double double double double HC_X HC_Y H_R2 SH_Y = = = = 1.78 * D / 2.; 3. * D; D * D / 4.; 1. * D; const double GAS_R = 8.31446261815324; const int const double const double COMP_COUNT = 2; COMP_CP[COMP_COUNT] = { 1006., 1006. }; COMP_M[COMP_COUNT] = { 28.98e-3, 28.98e-3 }; const const const const const const const double R0 = 1.65, R1 = 1.20, R2 = 0.166; double U0 = 0.0, U1 = 0.0, U2 = 0.0; double V0 = 114.4, V1 = 0.0, V2 = 0.0; double P0 = 1.589e5, P1 = 1.01325e5, P2 = 1.01325e5; double C0[COMP_COUNT] = { 0.,1. }; double C1[c OMP_COUNT] = { 0.,1. }; double C2[c OMP_COUNt ] = { 1.,0. }; double** double** double** double** double** double** double** c[COMP_COUNT]; r; u; v; p; E; e; double** ru; double** rv; 47

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А double** re; double** rc[COMP_COUNT]; double** double** double** double** ru_old; rv_old; re_old; rc_old[COMP_COUNT]; void draw(char* fName, int nx, int ny); double KR(double x); double _max_(double x, double y); const const const const int int int int lo[2] = { K_WENO, K_WENO }; hi[2] = { K_WENO + NX - 1, K_WENO + NY - 1 }; NX_ALL = NX + K_WENO + K_WENO; NY_ALL = NY + K_WENO + K_WENO; char str[1024]; double calc_gam(int i, int j) { double cp_mix = 0.; double M_ = 0.; for (int ic = 0; ic < COMP_COUNT; ic++) { M_ += c[ic][i][j] / COMP_M[ic]; cp_mix += c[ic][i][j] * COMP_CP[ic]; } double M_mix = 1. / M_; double cv_mix = cp_mix - GAS_R / M_mix; return cp_mix / cv_mix; } void init() { for (int ic = 0; ic < COMP_COUNT; ic++) c[ic] = new double* [NX_ALL]; r = new double* [NX_ALL]; u = new double* [NX_ALL]; v = new double* [NX_ALL]; p = new double* [NX_ALL]; E = new double* [NX_ALL]; e = new double* [NX_ALL]; ru = new rv = new re = new for (int double* double* double* ic = 0; [NX_ALL]; [NX_ALL]; [NX_ALL]; ic < COMP_COUNT; ic++) rc[ic] = new double* [NX_ALL]; ru_old = rv_old = re_old = for (int new double* [NX_ALL]; new double* [NX_ALL]; new double* [NX_ALL]; ic = 0; ic < COMP_COUNT; ic++) rc_old[ic] = new double* [NX_ALL]; for (int i = 0; i < NX_ALL; i++) { for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i] = new double[NY_ALL]; r[i] = new double[NY_ALL]; u[i] = new double[NY_ALL]; v[i] = new double[NY_ALL]; 48

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А p[i] = new double[NY_ALL]; E[i] = new double[NY_ALL]; e[i] = new double[NY_ALL]; ru[i] = new rv[i] = new re[i] = new for (int ic ru_old[i] = rv_old[i] = re_old[i] = for (int ic double[NY_ALL]; } double[NY_ALL]; double[NY_ALL]; double[NY_ALL]; = 0; ic < COMP_COUNT; ic++) rc[ic][i] = new double[NY_ALL]; new double[NY_ALL]; new double[NY_ALL]; new double[NY_ALL]; = 0; ic < COMP_COUNT; ic++) rc_old[ic][i] = new for (int i = lo[0]; i <= hi[0]; i++) { double x = XMIN + (i - lo[0]) * HX; for (int j = lo[1]; j <= hi[1]; j++) { double y = YMIN + (j - lo[1]) * HY; if (y < SH_Y) { r[i][j] = R0; //post shock u[i][j] = U0; v[i][j] = V0; p[i][j] = P0; for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i][j] = C0[ic]; } else { r[i][j] = R1; // pre shock u[i][j] = U1; v[i][j] = V1; p[i][j] = P1; for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i][j] = C1[ic]; } double gam = calc_gam(i, j); e[i][j] = p[i][j] / r[i][j] / (gam - 1.0); E[i][j] = e[i][j] + (u[i][j] * u[i][j] + v[i][j] * v[i][j]) / 2.0; ru[i][j] = r[i][j] * u[i][j]; rv[i][j] = r[i][j] * v[i][j]; re[i][j] = r[i][j] * E[i][j]; for (int ic = 0; ic < COMP_COUNT; ic++) rc[ic][i][j] = r[i][j] * c[ic][i][j]; } } } void bnd_cond() { 49

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А for (int j = lo[1]; j { for (int i = 0; { int i_ex int i_in <= hi[1]; j++)//no Y(left and right borders) i < K_WENO; i++) = lo[0] - i - 1; = lo[0] + i; r[i_ex][j] = r[i_in][j]; for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i_ex][j] = c[ic][i_in][j]; u[i_ex][j] v[i_ex][j] p[i_ex][j] e[i_ex][j] E[i_ex][j] = = = = = -u[i_in][j]; v[i_in][j]; p[i_in][j]; p[i_ex][j] / r[i_ex][j] / (calc_gam(i_in, j) - 1.0); e[i_ex][j] + (u[i_ex][j] * u[i_ex][j] + v[i_ex][j] * v[i_ex][j]) * 0.5; for (int ic r[i_ex][j] * c[ic][i_ex][j]; ru[i_ex][j] rv[i_ex][j] re[i_ex][j] = 0; ic < COMP_COUNT; ic++) rc[ic][i_ex][j] = =r[i_ex][j] =r[i_ex][j] =r[i_ex][j] * u[i_ex][j]; * v[i_ex][j]; * E[i_ex][j]; i_ex = hi[0] + i + 1; i_in = hi[0] - i; r[i_ex][j] = r[i_in][j]; for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i_ex][j] = c[ic][i_in][j]; u[i_ex][j] v[i_ex][j] p[i_ex][j] e[i_ex][j] E[i_ex][j] = = = = = -u[i_in][j]; v[i_in][j]; p[i_in][j]; p[i_ex][j] / r[i_ex][j] / (calc_gam(i_in, j) - 1.0); e[i_ex][j] + (u[i_ex][j] * u[i_ex][j] + v[i_ex][j] * v[i_ex][j]) * 0.5; for (int ic = 0; ic < COMP_COUNT; ic++) rc[ic][i_ex][j] = r[i_ex][j] * c[ic][i_ex][j]; ru[i_ex][j] =r[i_ex][j] * u[i_ex][j]; rv[i_ex][j] =r[i_ex][j] * v[i_ex][j]; re[i_ex][j] =r[i_ex][j] * E[i_ex][j]; } } for (int i = lo[0]; i <= hi[0]; i++)//no X(top and bottom borders) { for (int j = 0; j < K_WENO; j++) { int j_ex = lo[1] - j - 1; int j_in = lo[1] + j; //Реализовано постоянное втекание(на нижней границе всегда постоянное значение POSTshock) //и вытекание (на верхней гранце скорость газа не отражается) r[i][j_ex] = R0; for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i][j_ex] = C0[ic]; u[i][j_ex] = U0; v[i][j_ex] = V0; p[i][j_ex] = P0; e[i][j_ex] = p[i][j_ex] / r[i][j_ex] / (calc_gam(i, j_in) - 1.0); 50

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А E[i][j_ex] = e[i][j_ex] + (u[i][j_ex] * u[i][j_ex] + v[i][j_ex] * v[i][j_ex]) * 0.5; for (int ic = 0; ic < COMP_COUNT; ic++) rc[ic][i][j_ex] = r[i][j_ex] * c[ic][i][j_ex]; ru[i][j_ex] = r[i][j_ex] * u[i][j_ex]; rv[i][j_ex] = r[i][j_ex] * v[i][j_ex]; re[i][j_ex] = r[i][j_ex] * E[i][j_ex]; j_ex = hi[1] + j + 1; j_in = hi[1] - j; r[i][j_ex] = r[i][j_in]; for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i][j_ex] = c[ic][i][j_in]; u[i][j_ex] v[i][j_ex] p[i][j_ex] e[i][j_ex] E[i][j_ex] = = = = = u[i][j_in]; v[i][j_in]; p[i][j_in]; p[i][j_ex] / r[i][j_ex] / (calc_gam(i, j_in) - 1.0); e[i][j_ex] + (u[i][j_ex] * u[i][j_ex] + v[i][j_ex] * v[i][j_ex]) * 0.5; for (int ic = 0; ic < COMP_COUNT; ic++) rc[ic][i][j_ex] = r[i][j_ex] * c[ic][i][j_ex]; ru[i][j_ex] = r[i][j_ex] * u[i][j_ex]; rv[i][j_ex] = r[i][j_ex] * v[i][j_ex]; re[i][j_ex] = r[i][j_ex] * E[i][j_ex]; } } } bool trash(int x_point, int y_point, int x_size, int y_size, int i, int j) { bool res; ((i >= x_point) && (i <= x_point + x_size) && (j >= y_point) && (j <= y_point + y_size)) ? res = true : res = false; return res; } //Отступ от левой и правой границы int BackStep = 40; /* //каверна==================================================================== const int x_point1 = lo[0] + 65; const int x_point2 = hi[0]; const int y_point1 = 160; //размер объекта const int x_size1 = x_point2 - x_point1;//длина const int y_size1 = hi[1]-y_point1;//ширина //вторая сторона каверны const int x_point1_2 = lo[0] + 65; const int x_point2_2 = hi[0]; const int y_point1_2 = lo[1]; //размер объекта const int x_size1_2 = x_point2_2 - x_point1_2;//длина const int y_size1_2 = 120-y_point1_2;//ширина 51

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А //========================================================================== */ //ступенька=============================================================== const int x_point1 = lo[0] + 65; const int x_point2 = hi[0]; const int y_point1 = 120; //размер объекта const int x_size1 = x_point2 - х_ро^^;//длина const int y_size1 = hi[1]-y_point1;//ширина //======================================================================== //два куба================================================================ //точка расположения объекта(нижнего левого угла) //const int x_point1 = lo[0] + BackStep; //const int x_point2 = hi[0] - BackStep; //const int y_point1 = 120; //размер объекта //const int x_size1 = x_point2 - x_point1;//длина //const int y_size1 = x_point2 - x_point1;//ширина //точка расположения объекта(нижнего левого угла) //const int x_point1_2 = lo[0] + BackStep; //const int x_point2_2 = hi[0] - BackStep; //const int y_point1_2 = 145; //размер объекта //const int x_size1_2 = x_point2_2 - x_point1_2;//длина //const int y_size1_2 = x_point2_2 - x_point1_2;//ширина //========================================================================= void lr_bnd_cond_cube(int x_point, int y_point, int x_size, int y_size) { for (int j = y_point; j <= y_point + y_size; ]++)//по Y(left and right borders) { for (int i = 0; i < K_WENO; i++) { int i_in = x_point + x_size + i + 1;//for right border int i_ex = x_point + x_size - i; r[i_ex][j] = r[i_in][j]; for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i_ex][j] = c[ic][i_in][j]; u[i_ex][j] v[i_ex][j] p[i_ex][j] e[i_ex][j] E[i_ex][j] = = = = = -u[i_in][j]; v[i_in][j]; p[i_in][j]; p[i_ex][j] / r[i_ex][j] / (calc_gam(i_in, j) - 1.0); e[i_ex][j] + (u[i_ex][j] * u[i_ex][j] + v[i_ex][j] * v[i_ex][j]) * 0.5; for (int ic r[i_ex][j] * c[ic][i_ex][j]; ru[i_ex][j] rv[i_ex][j] re[i_ex][j] = 0; ic < COMP_COUNT; ic++) rc[ic][i_ex][j] = =r[i_ex][j] * u[i_ex][j]; =r[i_ex][j] * v[i_ex][j]; =r[i_ex][j] * E[i_ex][j]; i_in = x_point - i - 1;//for left border i_ex = x_point + i; 52

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А r[i_ex][j] = r[i_in][j]; for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i_ex][j] = c[ic][i_in][j]; u[i_ex][j] v[i_ex][j] p[i_ex][j] e[i_ex][j] E[i_ex][j] = = = = = -u[i_in][j]; v[i_in][j]; p[i_in][j]; p[i_ex][j] / r[i_ex][j] / (calc_gam(i_in, j) - 1.0); e[i_ex][j] + (u[i_ex][j] * u[i_ex][j] + v[i_ex][j] * v[i_ex][j]) * 0.5; for (int ic = 0; ic < COMP_COUNT; ic++) rc[ic][i_ex][j] = r[i_ex][j] * c[ic][i_ex][j]; ru[i_ex][j] =r[i_ex][j] * u[i_ex][j]; rv[i_ex][j] =r[i_ex][j] * v[i_ex][j]; re[i_ex][j] =r[i_ex][j] * E[i_ex][j]; } } } void du_bnd_cond_cube(int x_point, int y_point, int x_size, int y_size) { for (int i = x_point; i <= x_point + x_size; i++)//no X(top and bottom borders) { for (int j = 0; j < K_WENO; j++) { int j_in = y_point + y_size + j + 1;//for top border int j_ex = y_point + y_size - j; r[i][j_ex] = r[i][j_in]; for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i][j_ex] = c[ic][i][j_in]; u[i][j_ex] v[i][j_ex] p[i][j_ex] e[i][j_ex] E[i][j_ex] = = = = = u[i][j_in]; -v[i][j_in]; p[i][j_in]; p[i][j_ex] / r[i][j_ex] / (calc_gam(i, j_in) - 1.0); e[i][j_ex] + (u[i][j_ex] * u[i][j_ex] + v[i][j_ex] * v[i][j_ex]) * 0.5; for (int ic = 0; ic < COMP_COUNT; ic++) rc[ic][i][j_ex] = r[i][j_ex] * c[ic][i][j_ex]; ru[i][j_ex] = r[i][j_ex] * u[i][j_ex]; rv[i][j_ex] = r[i][j_ex] * v[i][j_ex]; re[i][j_ex] = r[i][j_ex] * E[i][j_ex]; j_in = y_point - j - 1;//for bottom border j_ex = y_point + j; r[i][j_ex] = r[i][j_in]; for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i][j_ex] = c[ic][i][j_in]; u[i][j_ex] v[i][j_ex] p[i][j_ex] e[i][j_ex] E[i][j_ex] = = = = = u[i][j_in]; -v[i][j_in]; p[i][j_in]; p[i][j_ex] / r[i][j_ex] / (calc_gam(i, j_in) - 1.0); e[i][j_ex] + (u[i][j_ex] * u[i][j_ex] + v[i][j_ex] * v[i][j_ex]) * 0.5; 53

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А for (int ic = 0; ic < COMP_COUNT; ic++) rc[ic][i][j_ex] = r[i][j_ex] * c[ic][i][j_ex]; ru[i][j_ex] =r[i][j_ex] * u[i][j_ex]; rv[i][j_ex] =r[i][j_ex] * v[i][j_ex]; re[i][j_ex] =r[i][j_ex] * E[i][j_ex]; } } } double double double double double double double rl, pl, ul, vl, cl[COMP_COUNT]; rr, pr, ur, vr, cr[COMP_COUNT]; r_[2 * K_WENO]; c_[ c OMP_COUNt ][2 * K_WENO]; p_[2 * K_WENo ]; u_[2 * K_WENo ]; v_[2 * K_WENo ]; void WENO(double* UU, double& Um, double& Up); void calc_flux_hllc( double rl, double ul, double vl, double pl, double cl[], double gaml, double rr, double ur, double vr, double pr, double cr[], double gamr, double& qu, double& qv, double& qe, double qc[]); void step() { double thx = TAU / HX; double thy = TAU / HY; double fu, fv, fe, fc[COMP_COUNT]; // X direction lr_bnd_cond_cube(x_point1, y_point1, x_size1, y_size1); //lr_bnd_cond_cube(x_point1_2, y_point1_2, x_size1_2, y_size1_2); for (int i = lo[0] - 1; i <= hi[0]; i++) { for (int j = lo[1]; j <= hi[1]; j++) { for (int k = -K_WENO + 1; k <= K_WENO; k++) { for (int ic = 0; ic < COMP_COUNT; ic++) c_[ic][k + K_WENO 1] = c[ic][i + k][j]; r_[k + K_WENO - 1] = r[i + k][j]; p_[k + K_WENO - 1] = p[i + k][j]; u_[k + K_WENO - 1] = u[i + k][j]; v_[k + K_WENO - 1] = v[i + k][j]; } for (int ic = 0; ic < COMP_COUNT; ic++) WENO(c_[ic], cl[ic], cr[ic]); WENO(r_, rl, rr); WENo(p_, pl, pr); WENO(u_, ul, ur); WENo(v_, vl, vr); double gaml = calc_gam(i, j); double gamr = calc_gam(i + 1, j); calc_flux_hllc(rl, ul, vl, pl, cl, gaml, rr, ur, vr, pr, cr, gamr, 54

ss 'Л44. * эх =+ [X + C][T]aj * aj. =+ [x + C][t ]aj * nj- =+ [I + P] [T]nu =+ [X + Р][Т][эт]эи (++ЭТ fINflOD-dWOD > эт '0 = эт хит) { 1Лчх { * [эт]эх jo j . 1Лцх * aj. =- [С][т]эи 1Лцх * aj. =- [ C ] [ t ] aj 1Лцх * nj- =- [P][T]nj * [эт]эх =- [P][t ][d t ]dj (++ЭТ fINHOD-dWOD > эт '0 = эт хит) 'ХшвЗ 1Лцх jox •(3J- 'эх 'nx 'лх 'juibS 'jd 'jd 'un 'ja 'jj fx3 'id fIn 'тл fxj )3II4-xnIJ--3IB3 P 'т)шв2- эхвэ = juibS axqnop f(C 'т)швЗ- эхвэ = хшв§ axqnop f(X + 1(ил 'хл '_A)ON3M f(jn fxn '_n)ON3M f(jd fxd f“d)0N3M f(jj fxu f—u)ON3M 1 ([эт]иэ '[эт]хэ f[:>T] 3)0N3M (++эт fINHOD dWOD > эт '0 = эт хит) ‘•[>\ + ?[>l + ?[>l + Р][т]л = [x - 0N3M_X + >1]_л P][T]n = [x - 0N3M->I + J|]-n P][l]d = [X - 0N3M- >1 + >|]“d jox { J[>l + Р][Т][эт]э = [x - 0N 3M >1 + >|][эт] э (++эт fINHOD dWOD > эт '0 = d t j |.ut) jox ?[>l + P][T]J = [X - ON3M->1 + >|]_J } (++>l ?ON3M >1 => >1 ‘ X + ON3M X- = >1 * ut) jox } (++C f [Х]ТЧ => P ‘X - [x]°x = p хит) jox (++t } f[0 ]xq => т f[0 ]ox = т хит) jox f(3 - X3ZTS- A f3- X9ZTS- x 'z - TXUTod- A 'z - TXUTod- x)aqno- puoo- puq- n p / / f(xazTS- A fxazTS- x 'ixuT od- A 'TXUTod- x)aqno- puoo- puq- np иотхээитр л // fxqx * fxqx * fxqx * эх =+ [C ][x + т]эи [C ][x + t ] a j nj-=+ [ P][X + T]n j { { ax =+ fxqx * [эт]эх =+ [P][X + t ] [ ot ] oj (++эт fINHOD-dWOD > эт '0 = эт хит) jox fxqx * эх =- [С][т]эи fxqx * лх =- [P][t]aj fxqx * nj- =- [P][T]ru fxqx * [эт]эх =- [Р][т][эт]эи (++эт fINHOD-dWOD > эт f0 = эт хит) jox f(ox 'эх fлх 'nx V КИННЖОКИсШ эинэжиойodn

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А } void fix_c() { for (int i = lo[0]; i <= hi[0]; i++) { for (int j = lo[1]; j <= hi[1]; j++) { for (int ic = 0; ic < COMP_COUNT; ic++) { if (c[ic][i][j] < 0.) c[ic][i][j] = 0.; if (c[ic][i][j] > 1.) c[ic][i][j] = 1.; rc[ic][i][j] = r[i][j] * c[ic][i][j]; } } } } int main(int argc, char* argv[]) { init(); double t = 0.0; int i_step = 0; sprintf(str, "ro.%010d", i_step); draw(str, NX, NY); while (t < TMAX) { t += TAU; i_step++; bnd_cond(); for (int i = 0; i < NX_ALL; i++) { for (int ic = 0; ic < COMP_COUNT; ic++) memcpy(rc_old[ic][i], rc[ic][i], sizeof(double) * (NY_ALL)); memcpy(ru_old[i], ru[i], sizeof(double) * (NY_ALL)); memcpy(rv_old[i], rv[i], sizeof(double) * (NY_ALL)); memcpy(re_old[i], re[i], sizeof(double) * (NY_ALL)); } step(); for (int i = lo[0]; i <= hi[0]; i++) { for (int j = lo[1]; j <= hi[1]; j++) { r[i][j] = 0.; for (int ic = 0; ic < COMP_COUNT; ic++) r[i][j] += rc[ic][i][j]; for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i][j] = rc[ic][i][j] / r[i][j]; u[i][j] v[i][j] E[i][j] e[i][j] = = ru[i][j] / r[i][j]; = rv[i][j] / r[i][j]; = re[i][j] / r[i][j]; E[i][j] - (u[i][j] * u[i][j] + v[i][j] * v[i][j]) * 0.5; p[i][j] = r[i][j] * e[i][j] * (calc_gam(i, j) - 1 .0); 56

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А } } fix_c(); bnd_cond(); step(); for (int i = lo[0]; i <= hi[0]; i++) { for (int j = lo[1]; j <= hi[1]; j++) { for (int ic = 0; ic < COMP_COUNT; ic++) rc[ic][i][j] += 3. * rc_old[ic][i][j]; ru[i][j] += 3. * ru_old[i][j]; rv[i][j] += 3. * rv_old[i][j]; re[i][j] += 3. * re_old[i][j]; for (int ic = 0; ic < COMP_COUNT; ic++) rc[ic][i][j] /= 4.; ru[i][j] /= 4.; rv[i][j] /= 4.; re[i][j] /= 4.; r[i][j] = 0.; for (int ic = 0; ic < COMP_COUNT; ic++) r[i][j] += rc[ic][i][j]; for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i][j] = rc[ic][i][j] / r[i][j]; u[i][j] v[i][j] E[i][j] e[i][j] = = ru[i][j] / r[i][j]; = rv[i][j] / r[i][j]; = re[i][j] / r[i][j]; E[i][j] - (u[i][j] * u[i][j] + v[i][j] * v[i][j]) * 0.5; p[i][j] = r[i][j] * e[i][j] * (calc_gam(i, j) - 1.0); } } fix_c(); bnd_cond(); step(); for (int i = lo[0]; i <= hi[0]; i++) { for (int j = lo[1]; j <= hi[1]; j++) { for (int ic = 0; ru[i][j] *= 2.; rv[i][j] *= 2.; re[i][j] *= 2.; for (int ic = 0; rc_old[ic][i][j]; ru[i][j] += ru_o rv[i][j] += rv_o re[i][j] += re_o for (int ic = 0; ru[i][j] /= 3.; rv[i][j] /= 3.; re[i][j] /= 3.; r[i][j] = 0. ; 57

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А for (int ic = 0; ic < COMP_COUNT; ic++) r[i][j] += rc[ic][i][j]; for (int ic = 0; ic < COMP_COUNT; ic++) c[ic][i][j] = rc[ic][i][j] / r[i][j]; u[i][j] v[i][j] E[i][j] e[i][j] = = ru[i][j] / r[i][j]; = rv[i][j] / r[i][j]; = re[i][j] / r[i][j]; E[i][j] - (u[i][j] * u[i][j] + v[i][j] * v[i][j]) * 0.5; p[i][j] = r[i][j] * e[i][j] * (calc_gam(i, j) - 1.0); } } fix_c(); if (i_step % STEP_OUT == 0) { sprintf(str, "ro.%010d", i_step); draw(str, hi[0] - lo[0] + 1, hi[1] - lo[1] + 1); } if (i_step % 100 == 0) printf("step = %d\t\tt = %e\n", i_step, t); } for (int i = 0; { for (int delete[] delete[] delete[] delete[] delete[] delete[] i < NX_ALL; i++) ic = 0; ic < COMP_COUNT; ic++) delete[] c[ic][i]; r[i]; u[i]; v[i]; p[i]; E[i]; e[i]; for (int delete[] delete[] delete[] ic = 0; ic < COMP_COUNT; ic++) delete[] rc[ic][i]; ru[i]; rv[i]; re[i]; for (int delete[] delete[] delete[] ic = 0; ic < COMP_COUNT; ic++) delete[] rc_old[ic][i]; ru_old[i]; rv_old[i]; re_old[i]; } for (int delete[] delete[] delete[] delete[] delete[] delete[] ic = 0; ic < COMP_COUNT; ic++) delete[] c[ic]; r; u; v; p; E; e; for (int delete[] delete[] delete[] ic = 0; ic < COMP_COUNT; ic++) delete[] rc[ic]; ru; rv; re; for (int ic = 0; ic < COMP_COUNT; ic++) delete[] rc_old[ic]; 58

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А delete[] ru_old; delete[] rv_old; delete[] re_old; return 0; } double KR(double x) { if (x >= 0.01 && x <= 0.03) return -0.01 * sin(M_PI * (x - 0.01) / 0.02); return 0.0; } inline double _MAX_(double x, double y) { if (x <= y) return y; return x; } void WENO(double* UU, double& Um, double& Up) { double BETA[3]; double ALPHA[3]; double eps = 1.0e-6; if ((UU[2] - UU[1]) * (UU[3] - UU[2]) < 0.0) Um = UU[2]; else { BETA[0] = (13. / 12.) * (UU[2] - 2 * UU[3] + UU[4]) * (UU[2] - 2 * UU[3] + UU[4]) +0.25 * (3 * UU[2] - 4 * Uu [3] + UU[4]) * (3 * UU[2] - 4 * UU[3] + UU[4]); BETA[1] = (13. / 12.) * (UU[1] - 2 * UU[2] + UU[3]) * (UU[1] - 2 * Uu [2] + UU[3]) +0.25 * (UU[1] - UU[3]) * (UU[1]- UU[3]); BETA[2] = (13. / 12.) * (UU[0] - 2 * UU[1] + UU[2]) * (UU[0] - 2 * UU[1] + UU[2]) + 0.25 * (UU[0] - 4 * UU[1] + 3 * UU[2]) * (UU[0] - 4 * UU[1] + 3 * UU[2]); ALPHA[0] = 0.3 / ((eps + BETA[0]) * (eps + BETA[0])); ALPHA[1] = 0.6 / ((eps + BETa [1]) * (eps + BETa [1])); ALPHa [2] = 0.1 / ((eps + BETa [2]) * (eps + BETa [2])); Um = (ALPHA[0] * (2 * UU[2] + 5 * UU[3] - UU[4]) + ALPHA[1] * (-UU[1] + 5 * UU[2] + 2 * UU[3]) + ALPHA[2] * (2 * UU[0] - 7 * UU[1] + 11 * UU[2])) / ((ALPHA[0] + ALPHA[1] + ALPHA[2]) * 6); } if ((UU[3] - UU[2]) * (UU[4] - UU[3]) < 0.0) Up = UU[3]; else { BETA[0] = (13. / 12.) * (UU[3] - 2 * UU[4] + UU[5]) * (UU[3] - 2 * UU[4] + UU[5]) +0.25 * (3 * UU[3] - 4 * Uu [4] + UU[5]) * (3 * UU[3] - 4 * UU[4] + UU[5]); BETA[1] = (13. / 12.) * (u U[2] - 2 * UU[3] + UU[4]) * (UU[2] - 2 * Uu [3] + UU[4]) +0.25 * (UU[2] - UU[4]) * (UU[2]- UU[4]); BETA[2] = (13. / 12.) * (UU[1] - 2 * UU[2] + UU[3]) * (UU[1] - 2 * UU[2] + UU[3]) + 0.25 * (UU[1] - 4 * UU[2] + 3 * UU[3]) * (UU[1] - 4 * UU[2] + 3 * UU[3]); ALPHA[0] = 0.1 / ((eps + BETA[0]) * (eps + BETA[0])); ALPHA[1] = 0.6 / ((eps + BETa [1]) * (eps + BETA[1])); ALPHA[2] = 0.3 / ((eps + BETa [2]) * (eps + BETA[2])); Up = (ALPHA[0] * (11 * UU[3] - 7 * UU[4] + 2 * UU[5]) + ALPHA[1] * (2 * UU[2] + 5 * UU[3] - UU[4]) + ALPHA[2] * (-UU[1] + 5 * UU[2] + 2 * UU[3])) / ((ALPHA[0] + ALPHA[1] + ALPHA[2]) * 6); } } 59

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А void draw(char* fName, int nx, int ny) { char htmlName[50]; strcpy(htmlName, fName); strcat(htmlName, ".vtk"); FILE* fp = fopen(htmlName, "w"); fprintf(fp, "# vtk DataFile Version 2.0\n"); fprintf(fp, "results\n"); fprintf(fp, "ASCII\n"); fprintf(fp, "DATASET UNSTRUCTURED_GRID\n"); int pCount = (nx + 1) * (ny + 1); fprintf(fp, "POINTS %d float\n", pCount); for (int j = 0; j <= ny; j++) { for (int i = 0; i <= nx; i++) { fprintf(fp, "%f %f %f \n", XMIN + HX * i, YMIN + HY * j, 0.0); } } fprintf(fp, "\n"); int cellsCount = nx * ny; fprintf(fp, "CELLS %d %d\n", cellsCount, 5 * cellsCount); for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) { fprintf(fp, "4 %d %d %d %d \n", j * (nx + 1) + i, j * (nx + 1) + i + 1, (j + 1) * (nx + 1) + i + 1, (j + 1) * (nx + 1) + i); } } fprintf(fp, "\n"); fprintf(fp, "CELL_TYPES %d\n", cellsCount); for (int i = 0; i < cellsCount; i++) fprintf(fp, "9\n"); fprintf(fp, "\n"); fprintf(fp, "CELL_DATA %d\n", cellsCount); fprintf(fp, "SCALARS Density float 1\nLOOKUP_TABLE default\n"); for (int i = lo[0]; i <= hi[0]; i++) { for (int j = lo[1]; j <= hi[1]; j++) { fprintf(fp, "%f ", r[i][j]); } fprintf(fp, "\n"); } fprintf(fp, "SCALARS Pressure float 1\nLOOKUP_TABLE default\n"); for (int i = lo[0]; i <= hi[0]; i++) { for (int j = lo[1]; j <= hi[1]; j++) { //will cut "trash" pressure in Paraview 60

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А (trash(x_point1, y_point1, x_size1, y_size1, i, j))? //|| trash(x_point1_2, y_point1_2, x_size1_2, y_size1_2, i, j) ) ? fprintf(fp, "%f ", 0) : fprintf(fp, "%f ", p[i][j]); if (p[i][j] != p[i][j]) { int zhrv = 0; } } fprintf(fp, "\n"); } fprintf(fp, "SCALARS Energy float 1\nLOOKUP_TABLE default\n"); for (int i = lo[0]; i <= hi[0]; i++) { for (int j = lo[1]; j <= hi[1]; j++) { fprintf(fp, "%f ", E[i][j]); } fprintf(fp, "\n"); } fprintf(fp, "SCALARS Mach_number float 1\nLOOKUP_TABLE default\n"); for (int i = lo[0]; i <= hi[0]; i++) { for (int j = lo[1]; j <= hi[1]; j++) { fprintf(fp, "%f ", sqrt((u[i][j] * u[i][j] + v[i][j] * v[i][j]) / (calc_gam(i, j) * p[i][j] / r[i][j]))); } fprintf(fp, "\n"); } fprintf(fp, "VECTORS Velosity float\n"); for (int i = lo[0]; i <= hi[0]; i++) { for (int j = lo[1]; j <= hi[1]; j++) { fprintf(fp, "%f %f %f ", u[i][j], v[i][j], 0.0); } fprintf(fp, "\n"); } for (int ic = 0; ic < COMP_COUNT; ic++) { fprintf(fp, "SCALARS C_%03d float 1\nLOOKUP_TABLE default\n", ic); for (int i = lo[0]; i <= hi[0]; i++) { for (int j = lo[1]; j <= hi[1]; j++) { fprintf(fp, "%f ", c[ic][i][j]); } fprintf(fp, "\n"); } } fclose(fp); printf("File '%s' saved...\n", htmlName); } 61

П р о д о л ж ен и е П Р И Л О Ж Е Н И Я А #define F_HLLC_U(UK, FK, SK, SS, PK, RK, VK) (((SS)*((SK)*(UK)-(FK)) + (SK)*( (PK)+(RK)*((SK)-(VK))*((SS)-(VK)) )) / ((SK)-(SS))) #define F_HLLC_V(UK, FK, SK, SS, PK, RK, VK) (((SS)*((SK)*(UK)-(FK))) / ((SK)-(SS))) #define F_HLLC_e (u K, FK, SK, SS, PK, RK, V k ) (((s s )*((s k )*(u k )-(f k )) + (SK)*( (PK)+(RK)*((SK)-(VK))*((SS)-(VK)) )*(SS)) / ((SK)-(SS))) void calc_flux_hllc( double rl, double ul, double vl, double pl, double cl[], double gaml, double rr, double ur, double vr, double pr, double cr[], double gamr, double& qu, double& qv, double& qe, double qc[]) { int i; double sl, sr, p_star, s_star, p_pvrs, ql, qr, tmp; double czl = sqrt(pl * gaml / rl); double czr = sqrt(pr * gamr / rr); double el = pl / (rl * (gaml - 1.)); double er = pr / (rr * (gamr - 1.)); double e_tot_l = el + 0.5 * (ul * ul + vl *vl); double e_tot_r = er + 0.5 * (ur * ur + vr *vr); p_pvrs = 0.5 * (pl + pr) - 0.5 * (ur - ul) * 0.25 * (rl + rr) * (czl + czr); p_star = (p_pvrs > 0.) ? p_pvrs : 0.; ql = (p_star <= pl) ? 1 : sqrt(1. + (gaml + 1.) * (p_star / pl - 1.) / (2. * gaml)); qr = (p_star <= pr) ? 1 : sqrt(1. + (gamr + 1.) * (p_star / pr - 1.) / (2. * gamr)); sl = ul - czl * ql; sr = ur + czr * qr; if (sl > sr) { tmp = sl; sl = sr; sr = tmp; } s_star s_star s_star s_star = pr - pl; += rl * ul * (sl - ul); -= rr * ur * (sr - ur); /= (rl * (sl - ul) - rr * (sr - ur)); if (s_star < sl) s_star = sl; if (s_star > sr) s_star = sr; if (!((sl <= s_star) && (s_star <= sr))) { printf("HLLC: inequaluty SL <= S* <= SR is FALSE.\n"); exit(1); } if (sl >= 0.) { qu = rl * ul * ul + pl; qv = rl * vl * ul; qe = (rl * e_tot_l + pl) * ul; 62

О кончание П Р И Л О Ж Е Н И Я А for (i = 0; i < COMP_COUNT; i++) { qc[i] = rl * ul * cl[i]; }} else if (sr <= 0.) { qu = rr * ur * ur + pr; qv = rr * vr * ur; qe = (rr * e_tot_r + pr) * ur; for (i = 0; i < COMP_COUNT; i++) { qc[i] = rr * ur * cr[i]; }} else { if (s_star >= 0) { qu = F_HLLC_U( /* UK, FK, SK, SS, PK, RK, VK */ rl * ul, rl * ul * ul + pl, sl, s star, pl, rl, ul ); qv = F_HLLC_V( /* UK, FK, SK, SS, PK, RK, VK */ rl * vl, rl * ul * vl , sl, s star, pl, rl, ul ); qe = F_HLLC_E( /* UK, FK, SK, SS, PK, RK, VK */ rl * e_tot_l, (rl * e_tot_l + pl) * ul, sl, s_star, pl, rl, ul ); for (i = 0; i < COMP_COUNT; i++) { qc[i] = F_HLLC_V( /* UK, FK, SK, SS, PK, RK, VK */ rl * cl[i], rl * cl[i] * ul, sl, s_star, pl, rl, u);}} else { qu = F_HLLC_U( /* UK, FK, SK, SS, PK, RK, VK */ rr * ur, rr * ur * ur + pr, sr, s_star, pr, rr, ur ); qv = F_HLLC_V( /* UK, FK, SK, SS, PK, RK, VK */ rr * vr, rr * ur * vr, sr, s_star, pr, rr, ur ); qe = F_HLLC_E( /* UK, FK, SK, SS, PK, RK, VK */ rr * e_tot_r, (rr * e_tot_r + pr) * ur, sr, s_star, pr, rr, ur ); for (i = 0; i < COMP_COUNT; i++) { qc[i] = F_HLLC_V( /* UK, FK, SK, SS, PK, RK, VK */ rr * cr[i], rr * cr[i] * ur, sr, s_star, pr, rr, ur ); }}}} 63

Заявление о самостоятельном характере выполнения выпускной квалификационной работы Я, Багапов Ариф Ренатович, студент 4 курса, направления подготовки Прикладная математика и информатика, заявляю, что в моей работе на тему «Моделирование акустических полей при обтекании тел потоком газа», представленной в Г осударственную экзаменационную комиссию для публичной защиты, не содержится элементов неправомерных заимствований. Все прямые заимствования из печатных и электронных источников, а также ранее защищенных письменных работ, кандидатских и докторских диссертаций имеют соответствующие ссылки. Я ознакомлен с действующим в Университете Положением о проверке работ обучающихся ФГБОУ ВО «МГУ им. Н.П. Огарева» на наличие заимствований, в соответствии с которым обнаружение неправомерных заимствований является основанием для отрицательного руководителя работы. Ф И О / Z Locujl 2 0 2 0 ? дат а подпись отзыва

Заведующему кафедрой прикладной математики, дифференциальных уравнений и теоретической механики, канд. физ.-мат. наук, доц., Р. В. Жалнину студента 4 курса очной формы обучения (на бесплатной основе) направления подготовки «Прикладная математика и информатика» факультета математики и информационных технологий Вагапова Арифа Ренатовича заявление Прошу разместить мою выпускную квалификационную работу на тему «Моделирование акустических полей при обтекании тел потоком газа» в электронной библиотечной системе университета в полном объеме. Ф И О // 0LH9HJL дат а 2 02 О? подпись

отзыв о выпускной квалификационной работе студента Вагапова Арифа Ренатовича обучающегося по направлению подготовки 01.03.02 «Прикладная математика и информатика» на тему «Моделирование акустических полей при обтекании тел потоком газа» Задача моделирования акустических полей при обтекании различных тел потоком газа является одной из актуальных задач современной прикладной математики. Наличие эффективных инструментов математического моделирования и оценки акустических характеристик может позволить конструировать изделия с минимальным акустическим воздействием на окружающие объекты. Это делает актуальным развитие методов моделирования акустических полей. Работа Вагапова А.Р. посвящена построению вычислительного алгоритма для моделирования акустических полей при обтекании тел потоком газа. В ходе работы автором самостоятельно были решены следующие задачи: - изучена литература и выполнен обзор современного состояния вопроса; - реализован алгоритм расчета акустических характеристик; - выполнено моделирование акустических полей при обтекании различных конфигураций тел. Во время выполнения работы Вагапов А.Р. продемонстрировал самостоятельность в решении поставленных перед ним задач, проявил усердие и устремленность к достижению поставленных целей. Бакалаврская работа Вагапова Арифа Ренатовича соответствует заданию и всем требованиям, предъявляемым к работам подобного рода. Считаю, что работа заслуживает оценки «отлично», а ее автор заслуживает присвоения степени бакалавра по направлению подготовки «Прикладная математика и информатика» Научный руководитель к.ф.-м.н., доцент Р.В. Жалнин

ОТЧЕТ о результатах проверки работы обучающегося на наличие заимствований Ф.И.О. автора работы Багапов Ариф Ренатович Тема работы Моделирование акустических полей при обтекании тел потоком газа Руководитель Жалнин Руслан Викторович Представленная работа прошла проверку на наличие заимствований в системе «Антиплагиат. ВУЗ» Результаты автоматической проверки: оригинальность 82.08 % цитирования 3.86 % заимствования 14.06 % Результаты анализа полного отчета на наличие заимствований: правомерные заимствования: 14.06 %____________________________ корректные цитирования: 3.86 %________________________________ неправомерные заимствования: нет______________________________ признаки обхода системы: нет__________________________________ Общее заключение об итоговой оригинальности работы и возможности ее допуска к защите: Студент допускается к предварительной защите и защите ВКР в ГЭК. Руководитель работы к.ф.-м.н., доцент Р. В. Жалнин

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв