Численное моделирование динамики волновых систем на основе явной схемы Мак-Кормака Текст научной статьи по специальности «Физика»

УДК 539.37

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ВОЛНОВЫХ СИСТЕМ НА ОСНОВЕ ЯВНОЙ СХЕМЫ МАК-КОРМАКА

Д.А. ГУБАЙДУЛЛИН, Д.А. ТУКМАКОВ ИММ КазНЦ РАН

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

Ключевые слова: акустический резонатор, уравнения Навье-Стокса, явная схема Мак-Кормака, нелинейные и разрывные колебания, затухание, добротность.

Введение

Теоретическое исследование колебаний столба вязкого газа, связанное с решением системы уравнений Навье-Стокса [1-5], представляет собой сложную задачу в силу того, что близость результатов численного и физического эксперимента, особенно при моделировании динамики систем с малым затуханием, определяется диссипативными свойствами численной схемы. В данной статье сопоставляются результаты численного решения явным конечно-разностным методом Мак-Кормака [3] ряда задач акустики газа с известными из литературы аналитическими и экспериментальными результатами. При численном моделировании важно выяснить, насколько точно метод описывает снижение резонансной частоты вследствие диссипации, вызванной вязкостью [3], включающей в себя как численную [6], так и физическую составляющие. В этом случае можно также сопоставить параметры добротности физической системы и ее математической модели и решить вопрос о границах применимости метода моделирования при описании колебательной системы.

Численный метод

Движение газа описывается системой уравнений Навье-Стокса в консервативной форме в осесимметричной подвижной системе координат [2- 5]:

% + О; (1)

г -Г

р ри ру е

4 = 1' 1'1'1

© Д.А. Губайдуллин, Д.А. Тукмаков Проблемы энергетики, 2012, № 5-6

Е =

5гР + 5 хРи + 5 у РУ 1

5Ри + 5х (ри2 + Р-Тхх ) + 5у ( - Тху )

1

5РУ + 5х (риу-тху ) + 5у ( + Р-Туу )

1

р =

5е + 5х ( + хх ху + Ох ) + 5у (( + уу ху + Оу )

1 .

Пг Р + ПхРи + П yРV 1

Пг Ри + Пх (ри 2 + Р -тхх) + Пу (

( -Т ху )

1

Пг РУ + Пх ( - Т ху) + П у ( + Р - Т уу) 1

Пге + Пх (( + Р -Тхх)и-ТxyV + °х) + Пу ((е + Р-Туу )-тхуи + °у)

1

о =

РУ / у 1

(-РиУ + Тху ) у 1

(-Ру2 + Т уу ) у 1

( + Р -туу )- к -ду + и Т ху ^ / у

ди ду

1

• Т1ГИ = Тл

Тху = J; ^ух ^ху

, ди 2 , ,тт. , ду 2 ,тт. . т

Т хх = 2Ц-—<Ми);Т уу = 2Ц—<Ми); Р = (у-1)р/ ;

дх 3

ду 3'

ди ду V 2 2

= -+ — + - ; I = е-0,5р(и2 + V2);

дх ду у

1 =

5х 5 у 5г Пх П у Пг 0 0 1

; 5г =-хг5х - уг5у; Пг = -хгПх - уПу •

Здесь р, и, V, е, I, Т, Р — плотность, декартовы составляющие скорости, полная и внутренняя энергия, температура и давление газа; у, ц — постоянная адиабаты и

© Проблемы энергетики, 2012, № 5-6

динамическая вязкость; т^, туу, тху — составляющие тензора вязких напряжений; 1— якобиан преобразования координат при переходе от физической области в переменных (х, у, ?) к расчетной области в переменных (Е, п , ) Система (1) решалась при помощи явного метода Мак-Кормака второго порядка точности в обобщенных подвижных координатах [2,3, 4]: Е = Е(х, у, ) п =п(х, у.

Коррекция решения, найденного методом Мак-Кормака, проводилась с помощью схемы нелинейной коррекции, описанной в работе [7] в применении к векторам газодинамических функций в примитивных переменных и=(р, и, V, е) . Процедура выполнялась по всем строкам (вдоль Е), а затем по всем столбцам (вдоль п) в расчетной области [4].

Результаты расчетов

Решение задачи о распаде разрыва. Цилиндрическая труба, закрытая с обоих концов разделена на два отсека диафрагмой. В обоих отсеках расположен одинаковый газ, имеющий различные давления. Газ, находящийся в левом отсеке сжат до значительно большего давления, чем газ, находящийся в правом отсеке. В некоторый момент времени диафрагма разрушается, и газ из левого отсека устремляется в правый. Разрыв давлений, имевший место до разрушения диафрагмы, распространяется в виде ударной волны вправо, увлекая за собой спутный поток газа. Аналитический расчет числа Маха М, скорости распространения ударной волны в и скорости спутного потока V при заданном перепаде давления выполняется по следующим формулам [8]:

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

В начальный момент времени во внутренних узлах расчетной области определялись температура и плотность неподвижного газа: Т= 293К, р = 1,3 кг/м . Начальное давление слева от поверхности разрыва равно произведению интенсивности разрыва на давление в правой части трубы, где оно составляет величину р1=0,714 в безразмерных единицах. Конечно-разностная сетка содержит Щ =200 узлов в осевом и Nk =4 узла в радиальном направлении.

На рис. (1, а, б) показаны волны давления и скорости спутного потока газа, полученные при численном и аналитическом решении задачи о распаде разрыва интенсивности р2/р1=1,2. Сплошной линией нанесены кривые давления и скорости, полученные численным расчетом, пунктирной линией обозначены кривые давления и скорости рассчитанные аналитически.

Расчеты показывают, что численное и аналитическое решения для волн давления и скорости спутного потока близки и изменяются одинаковым образом при увеличении интенсивности разрыва. Начальный скачок давления распадается на волну сжатия, распространяющуюся вправо, и волну разрежения, идущую в область повышенного давления — влево. В средней части волны образуется «полка», в пределах которой давление равно среднему арифметическому в волнах сжатия и разрежения. На рис. 1, б приводится сравнение численного и аналитического решения для скорости

(2)

спутного потока. Величина скорости спутного потока V, полученная численно, близка к найденной аналитически по формулам (2). В численной модели скорость ударной волны для перепада давления р2/р1=1,2 составляет 9ч=1,0744 с, аналитически рассчитанная скорость ударной волны 9а=1,0823с. Скорость спутного потока газа за фронтом ударной волны, рассчитанная аналитически и численно составляет Уа =21,96 м/с и Уч =21,8 м/с. Скорость спутного потока, полученная численно была измерена по средней точке фронта волны сжатия. Скорость волны разрежения в численных расчетах 0,95с. На границах фронта волны наблюдается незначительная осцилляция (рис. 1, б).

Р 0,860,840,820,800,780,76 0,74 0,72 0,70

0,0 0,2 0,4 0,6 0,8

1, Х,м

0,0 0,2 0,4 0,6

а б

Рис.1. Решение задачи о распаде разрыва в момент времени 1=0,22: а - давление, б - скорость спутного потока

Х,м

Решение задачи о формировании ударной волны на поршне. Поршень, расположенный на левом конце трубы, мгновенно разгоняется до скорости и0, после чего движется равномерно. Поршень толкает перед собой слой газа, в котором формируется ударная волна, бегущая по неподвижному невозмущенному газу с некоторой скоростью 9, оставляющая за собой возмущенный газ, выведенный из состояния покоя и движущийся со скоростью спутного потока У. Аналитический расчет скорости распространения ударной волны 0 и скорости спутного потока У проводят по следующим формулам [8]:

( \

9 =

У-1 + 7 + 1 Рг .с ' 2У 2У Р1 '

У = -

у +1

п 1 +Т + 1 Р2

1

2У 2У Р1

у-1 +_у +1 Р2

с.

(3)

2У 2У Р1 ^

Здесь р1 и р2 - давление невозмущенного и возмущенного газа. Начальная температура и плотность невозмущенного воздуха: Т = 293К, р = 1,3 кг/м ; начальное давление в безразмерных единицах р\ = 0,714; скорость поршня 17м/с. Задача решалась в системе координат с неподвижным поршнем, на поверхности которого составляющие скорости полагались равными нулю, а для остальных газодинамических переменных задавались однородные граничные условия второго рода. В начальный момент времени во внутренних узлах расчетной области задавались температура, плотность и скорость газа и(Х = 0) = -и0 , где и0 - скорость поршня в системе координат с движущимся поршнем. Конечно-разностная сетка содержала Щ = 200 узлов в осевом и Nk = 4 узла в радиальном направлении. На рис. 2 изображен процесс формирования волн скорости газа для ряда безразмерных моментов времени 0,112; 0,337; 0,563;

0,789. Там же последовательно изображен профиль развивающейся ударной волны. На рис. 3 сопоставлено численное решение и аналитическое решение, полученное по формуле (3); сравнение проведено для скорости спутного потока в момент времени t = 0,45. Аналитическое решение обозначено пунктирной линией, численное решение -непрерывной. Количественно и качественно они близки.

V,м/с

16

V, м/с

12

0,0

0,2

0,4

0,6

1 Х,м

0,0

0,2

0,4

0,6

0,8

Хм

Рис.2. Скорость спутного потока при формировании ударной волны при скорости поршня и0=17 м/с

Рис.3. Скорость спутного потока в момент времени г=0.45 при скорости поршня и0=17 м/с

Колебания вязкого газа в закрытой трубе. Численное моделирование продольных колебаний воздушного столба в закрытой трубе, генерируемых перемещающимся по гармоническому закону с амплитудой а = 0,0003 м поршнем, проводилось на частоте первого линейного резонанса f1 = с/(2Ь) и в его окрестности. Длина трубы, имеющей диаметр d = 0,0375 м, составляла Ь = 1,06 м. Задача решалась в осесимметричной постановке. На оси трубы ставились граничные условия симметрии для составляющих скорости и однородные граничные условия второго рода для остальных газодинамических функций, а на твердых поверхностях, в том числе и на поверхности движущегося поршня, для составляющих скорости задавались условия прилипания; для плотности, давления, энергии, температуры - однородные граничные условия второго рода. На неподвижных твердых поверхностях для составляющих скорости задавались однородные граничные условия первого рода, для остальных газодинамических функций - однородные граничные условия второго рода. На поверхности поршня задавался гармонический закон изменения осевой составляющей скорости, радиальная составляющая скорости приравнивалась нулю.

В начальный момент времени во внутренних узлах расчетной области задавались температура и плотность неподвижного газа: Т = 293К, р =1,3 кг/м . Конечно-разностная сетка содержала N = 200 узлов в осевом и N = 90 узлов в радиальном направлении.

Проводилось сравнение результатов, полученных численно с результатами физического эксперимента [9,10], в котором определялся безразмерный размах колебаний давления Ар = (р2-р1)/р0 на различных частотах. Здесь р2, р\, р0 -максимальное, минимальное и среднее значения давления за период. На рис. 4 изображены кривые зависимости давления от времени на закрытом конце трубы для колебаний на частотах, лежащих в окрестности первой собственной частоты ^: ^0,9-f= 1,09-/[. На частотах f= 0,9-f= 1,09-^ наблюдаются биения [11], возникающие в связи с генерацией в системе собственной частоты колебаний и затухающие с течением времени. Частота биений равна разности между частотой колебаний поршня и собственной частотой акустической системы, а время их существования определяется ее добротностью. © Проблемы энергетики, 2012, № 5-6

8

4

0

0

Определив частоту биений и зная частоту колебаний поршня, можно найти резонансную частоту системы [11]. Для частот внешнего воздействия, меньших резонансной, частота резонанса определяется как / = 1/Т +/ , где Т - период биения, / частота первого линейного резонанса в моделируемой колебательной системе , / частота воздействия внешней силы. Для частот внешнего воздействия больших резонансной, частота резонанса определяется как / = /-1/Т. Оценим величину первой собственной частоты акустической системы с затуханием, зная параметры биений генерируемых в системе на дорезонансной и зарезонансной частотах. При колебаниях поршня с частотой 0,9/1, где / = 161,7 Гц, генерируются биения с периодом Т = 0,07 с, следовательно, /р = 158,18 Гц. В зарезонансной области частота колебаний поршня /= 1,09/1, период биения Т = 0,055 с. Резонансная частота /р = /-1/7=158,3 Гц. Таким образом, первая собственная частота акустической системы с диссипацией лежит в интервале 158,18 Гц </р < 158,3 Гц.

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

[12]

5 = ^®12 -ю р2 = 2я/1^1 - /р2/ /2 и добротности б = ю р / 28=/р /^2^/2 - /р2 ^

в физической и численной моделях, где Ю1 - циклическая частота первого линейного резонанса в идеальной системе, юр - в системе с затуханием. На основе анализа резонансных кривых в физическом эксперименте 8э = 0,11ю!; бэ = 4,5. В численном эксперименте наблюдается более высокое затухание и меньшая, чем в физическом эксперименте, добротность: 8ч = 0,2ш1; бч = 2,5. Резонансные кривые, полученные в расчетах и в эксперименте, близки вблизи резонанса и справа от него (рис. 5). В дорезонансной области резонансная кривая, полученная численно, лежит выше и левее экспериментальной. Характер резонансной кривой в численном эксперименте говорит о меньшей добротности, свойственной численной модели колебательной системы по сравнению с исходной физической системой.

р

0,7220,7200,7180,7160,7140,712 0,710

0,00 0,02

0,04

ДР-102 2

1,5

1

0,5 0

0,06 0,08 0,10 С с

0,9 1,0 1,05

ю/а>1

Рис. 4. Зависимость давления от времени при Рис. 5. Зависимость относительного размаха колебаниях поршня с частотами пунктирной колебаний давления газа от частоты при линией - /= 0,89/, сплошной линией - /= амплитуде колебаний поршня 0,0003 м: □ -1,09/1. давление, полученное в физическом

эксперименте, давление, полученное численно, изображено непрерывной линией. АР = (Рмакс-

Рмин)/Р0.

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

Работа выполнена по ГК №14.740.11.0351, а также при поддержке РФФИ (грант № 09-01-97013).

Summary

Comparison of wave properties of acoustic systems and their numerical models in problems of ideal and viscous gas is executed. For problems about disintegration of rupture and shock wave formation on the piston comparison of speeds of front of a wave of compression and speed of an accompanying stream of gas has been spent at the set pressure difference, process of formation of a shock wave numerical and analytically is described. The numerical model of oscillate motions of a column of viscous gas in a vicinity of first own frequency is investigated. On the continuance beat frequency arising in a vicinity of the first linear resonance of acoustic system with attenuation, its own frequency of fluctuations is defined. Comparison of factors of attenuation and good quality of physical oscillatory system and its numerical model is spent.

Keywords: the acoustic resonator, the equations of nave-stocks, the obvious scheme the mccormack, nonlinear and explosive fluctuations.

Литература

1. Ilgamov M. A., Zaripov R. G., Galiullin R. G., Repin V. B. Nonlinear oscillations of a gas in a tube // Appl. Mech. Rev. 1996. Vol. 49. No3. P.137-154.

2. Тукмаков А.Л. Зависимость механизма дрейфа твердой частицы в нелинейном волновом поле от ее постоянной времени и длительности прохождения волновых фронтов// ПМТФ. 2011. №4. С.105-106.

3. Флетчер К. Вычислительные методы в динамике жидкостей. Т.2. М.: Мир. 1991. 551 с.

4. Steger J. L. Implicit Finite-Difference Simulation of Flow about Arbitrary Two-Dimensional Geometries // AIAA J. 1978. Vol. 16, No 7. P. 679-686.

5. Ковеня В.М., Тарнавский Г.А., Черный С.Г. Применение метода расщепления в задачах аэродинамики. Новосибирск: Наука. Сибир. Отд-ние, 1990. 247 с.

6. Шокин Ю.И., Яненко Н.Н. Метод дифференциального приближения. Применение к газовой динамике. Новосибирск: Наука. Сиб. отд-ние, 1985. 364 с.

7. Жмакин А.И., Фурсенко А.А. Об одной монотонной разностной схеме сквозного счета// ЖВМ и МФ. 1980. Т.20. № 4. С.1021-1031.

8. Лойцянский Л.Г. Механика жидкости и газа. М.: Физматлит,1959. 784 с.

9. Губайдуллин Д.А., Зарипов Р.Г., Ткаченко Л.А. Нелинейные колебания аэрозоля вблизи резонанса в закрытой трубе в безударно волновом режиме.// Актуальные проблемы механики сплошной среды. К 20-летию ИММ КазНЦ РАН. Т.1-Казань: Фолиант,2011. 224 с.

10. Ткаченко Л.А. Зарипов Р.Г. Особенности нелинейных колебаний аэрозоля в закрытой трубе в безударно-волновом режиме // Вестник Нижегородского университета им. Н.И. Лобачевского. Механика жидкости и газа. 2011. № 4 (3). C. 1171-1173.

11. Горелик Г.С. Колебания и волны. М.:Физ. - мат. ГИЗ, 1959. 572 с.

12. Вахитов Я. Ш. Теоретические основы электроакустики и электроакустической аппаратуры. М.: Искусство, 1982. 415.

Поступила в редакцию

26 января 2012 г.

Губайдуллин Дамир Анварович - д-р физ.-мат. наук, член-корр. РАН, директор Института механики и машиностроения Казанского научного центра Российской академии наук (ИММ КазНЦ РАН). Тел.: 8 (843) 236-52-89. E-mail: gubajdullin@mail.knc.ru.

Тукмаков Дмитрий Алексеевич - аспирант Института механики и машиностроения Казанского научного центра Российской академии наук (ИММ КазНЦ РАН). Тел.: 8 (843) 272-44-92. E-mail: tukmndn@yandex.ru