CC BY
71
Том XXXVII
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
2 006
№ 4
УДК 534.21
532.525.011.5
О ПРИМЕНЕНИИ СХЕМЫ МАК-КОРМАКА ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ НЕЛИНЕЙНЫХ ЯВЛЕНИЙ ПРИ РАСПРОСТРАНЕНИИ ПЛОСКОЙ ЗВУКОВОЙ ВОЛНЫ В КАНАЛАХ
А. В. ЧИРИХИН
Определены условия применимости явной схемы Мак-Кормака для численного моделирования распространения плоской звуковой волны на основе уравнений Эйлера и Навье — Стокса. Проиллюстрированы характерные этапы трансформации первоначально гармонического колебания в последовательность слабых ударных волн в канале постоянного сечения и определена граница влияния вязкости и теплопроводности реального газа на условия формирования газодинамического разрыва. На примере канала типа сопла Лаваля при отсутствии течения показано влияние его формы. Для случая распространения звуковой волны по течению в сопле Лаваля продемонстрированы сохранение потока ее средней энергии при сохранении формы сигнала, близкой к гармонической, и переход сигнала в разрывное состояние.
Функционирование аэродинамических труб, а также различных технологических устройств, содержащих в конструкции компрессоры, вентиляторы или другие генераторы акустических колебаний, может сопровождаться взаимодействием таких колебаний и физических явлений в потоке, полезным или отрицательным в зависимости от конкретных практических целей. Так, в работе [1] на основе численного моделирования показано, что звуковые гармонические колебания
с амплитудой пульсации температуры ~ 1% могут приводить к интенсификации спонтанного фазового перехода в условиях типа камеры Вильсона и значительному увеличению количества капель конденсата. Возможность технологического применения подобного явления, если оно реализуется при течении конденсирующейся среды в сопле Лаваля, делает актуальным проведение соответствующего численного анализа. Очевидно, что при исследовании взаимодействия акустических пульсаций параметров потока и процесса его спонтанной конденсации приходится применять полные уравнения движения и единый численный алгоритм в отличие от традиционного подхода к задаче акустики движущейся среды [2, 3] и ее конечноразностного решения [4]. При этом практический интерес представляет оценка пригодности для численного моделирования распространения звуковой волны алгоритма, основанного на явной схеме Мак-Кормака и использованного в [5, 6] для расчета скачков конденсации. Сходная задача применительно к методу Годунова была ранее рассмотрена в [7]. В свою очередь, в качестве характерного газодинамического явления, удобного для тестирования расчетов, можно рассмотреть переход гармонической звуковой волны в последовательность слабых разрывов, для времени формирования которого
в [8] имеется аналитическое решение.
При распространении звуковой волны по потоку в сопле Лаваля практический интерес представляет вопрос о взаимодействии звуковой волны с основным течением и, в частности, о степени сохранения ее энергии.
1. Постановка задачи. Рассмотрим одномерное нестационарное движение вязкого теплопроводного газа в канале с контуром 5 = 5(х). В данном случае можно воспользоваться упрощенным вариантом системы уравнений (20), (22), (32) из монографии [9], которая в дивергентной форме будет иметь вид:
dU dV
dt dx
т0 тг _ с тл0
= G,
U = 5U0, V = 5V0, G = 5G ,
U0 =(р, pu, ре), V0 = (pu, P + pu2, puh),
Gо ! 0 p dln 5 , 4ц d u 4ц
dx 3Re dx 3Re 1
d2u
dx2
du
dx
цд2Т
(к-1) PrRe dxz
(1.1)
е = h - P/p, h =----T + 0.5u , P = p/к, p = pT,
к-1
Re = (p0c0L)/ Pr = (ц0ср c0 =(KRT0 )05.
Здесь t — время; x — координата; u — скорость; p — плотность среды; p — давление; T — температура; к — показатель адиабаты; h — полная удельная энергия среды; е — ее удельная внутренняя энергия; ц, п — коэффициенты вязкости и теплопроводности; Re, Pr — числа Рейнольдса и Прандтля соответственно; с0 — скорость звука; L — характерный линейный масштаб; Ср — удельная теплоемкость при р = const; R — газовая постоянная.
Уравнения движения представлены в безразмерном виде, при этом в качестве масштабов термодинамических параметров использованы их величины в исходном невозмущенном состоянии (или в состоянии изоэнтропического торможения) Р0, Р0, T), П0, Ц0; масштаб скорости С0, энергии С02 , времени Т = L/C0.
Для интегрирования системы (1.1) применялась явная схема Мак-Кормака совместно с простейшим вариантом монотонизатора [10]. При этом линейная координата х разбивается с постоянным шагом Ax на N + 1 узлов i = 0, ..., N. В узловых точках формируются начальные значения функций Ui, Vi. Для внутренних точек i = 1,., N - 1 решение при шаге по времени At получается последовательным применением соотношений:
U13 = Ui -At[(V+1 -V)/Ax-Gi], U23 = 0.5{(u13 +Ui)-At[(Fi -Vi-1 Г/
'Ax -G,
13
(1.2)
U1 =U23 + Ax(Ui+i -2Ui +U-1)/3.
Здесь степенью 1/3 отмечен результат работы предиктора, степенью 2/3 — результат работы корректора и единицей — окончательное решение после работы монотонизатора. Во всех расчетах шаг интегрирования по времени определялся по соотношению
At = 0.7 шт [Ах/(|и| + с0)г].
Для сохранения второго порядка аппроксимации при вычислении конечных разностей в узле / = N - 1 использовалась параболическая экстраполяция по внутренним точкам. При
/''"’О
вычислении производных от скорости и температуры в правых частях & использовались безразностные трехточечные формулы численного дифференцирования из [11].
В нулевом узле граничные условия формировались следующим образом. При расчете стационарного течения в сопле Лаваля на шаге п + 1 использовались условия сохранения
начальных значений полной энергии Но и энтропии ^ и равенство расхода до при / = 0 расходу в критическом сечении сопла д*, рассчитанному по параметрам на предыдущем временном слое п [12]:
(лр-к)п+1 =(PоP-к)t=0, (Н0)п+1 = (Н0)^ (д0)п+1 = (д*)п. (1.3)
Сравнение экспериментальных данных по скачкам конденсации при трансзвуковом течении влажного воздуха в крупномасштабной трубе Т-128 и результатов расчета в [5, 6] показало адекватность численного моделирования таких течений на основе алгоритма (1.2) совместно с условиями (1.3).
При моделировании распространения звуковых волн предполагалось, что в начальный момент времени (или с некоторого момента времени после установления стационарного решения)
в нулевом узле безразмерные статические параметры и скорость начинают изменяться согласно соотношениям для гармонического звукового колебания [8]:
р = р1 + Ар1 8т(2л/?), р = (р)1/к, и = и1 + Аи1 8т(2л/), Аи1 = —Р1(р1р1)-1/2 . (1.4)
к
Здесь /=/0т — безразмерная частота, тождественно равная числу Струхаля; Ар1 и Аи1 — безразмерные амплитуды давления и скорости; параметры и1 = 0, р1 = р1 = 1 для исходного состояния покоя или принимают соответствующие начальные значения при наличии течения. Все представленные ниже расчеты выполнены для газа с показателем адиабаты к = 1.4.
2. Некоторые примеры численного моделирования распространения звуковой волны в каналах. В работе [7] при обсуждении особенностей численного исследования нестационарных течений на основе метода Годунова отмечалось, что существенным моментом для режимов
с гармоническим изменением граничных условий является соотношение между шагом разностной сетки Ах и длиной волны такого процесса X = С0 / /0, которым определяется величина дисперсионных эффектов. В частности, если Ах недостаточно мало по сравнению с X, возможно проявление затухания возмущения по мере его распространения в цилиндрическом канале. Для устранения затухания в рассмотренной численной задаче авторы [7] выбирали шаг разностной сетки
с учетом условия Ах < (0.02 0.025)Х, которое было получено на основе специально
проведенных расчетов. Отмеченное обстоятельство обусловливает необходимость проведения соответствующего анализа и применительно к алгоритму (1.2). В связи с этим предварительно рассмотрим, как соотношение между Ах и X влияет на результаты численного моделирования на основе алгоритма (1.2) распространения гармонической звуковой волны в канале постоянного сечения, заполненного неподвижным невязким и нетеплопроводным газом. В данном случае правые части О0 в (1.1) будут равны нулю, а сама система сведется к уравнениям Эйлера. При этом в качестве характерного линейного масштаба естественно принять длину волны (Ь = X), что приведет к соответствию масштаба времени т периоду колебаний и тождеству / = 1 в граничных условиях (1.4).
0.8 ---------------1-------------1-------------1-------------1------X------1-------------1---------- Т- " "— 0
-3 -2.5 -2 -1.5 ^
Рис. 2
Так, в верхней части рис. 1 сплошной линией 1 показаны два периода гармонического изменения по времени статического давления р во входном сечении канала с безразмерной амплитудой Ар1 = 10-2. Ниже сплошной линией 2 показан результат расчета распространения такого сигнала при = Ах/Х = 0.025 в канале постоянного сечения, протяженность которого Х1 измеряется в длинах X. Здесь и далее на соответствующих рисунках по оси ординат отложена приведенная амплитуда сигнала 81 = (р - р1)/Ар1.
Расчет наглядно демонстрирует наличие сильного влияния дисперсии, приводящего на расстоянии Х1 ~ 60 к двукратному уменьшению амплитуды сигнала. Сплошной линией 3 вверху показано изменение сигнала по времени при Х1 = 60, а штриховой линией — гармонический сигнал той же амплитуды. Расчеты показывают, что в рассмотренных условиях частота сигнала сохраняется, а его форма близка к гармонической.
Если шаг разностной сетки уменьшить до значения = 0.003, то амплитуда сигнала будет воспроизведена с точностью ~ 1%, а при S = 0.001 получим воспроизведение амплитуды на расстоянии Х1 = 60 с точностью до 0.1%. Результат соответствующего расчета показан линией 4, а на характерном изменении формы сигнала остановимся позднее.
На рис. 2 в двух масштабах по ординате (левая ось — диапазон 0.8 1, кривые 1; правая ось —
диапазон 0 1, кривые 2) представлены приведенные амплитуды сигналов 81 при х1 = 30
О 10 20 30 40 50 X;
3 2 1
-1
37 38 39 40
-1
96
97
98
х-
Рис. 3
(штриховые линии) и Х1 = 60 (сплошные линии) при Др1 = 10 3 в зависимости от ^S, которые
пригодны для прогноза качества конкретных расчетов и оптимального выбора шага разностной сетки при других значениях протяженности канала. Для этого следует воспользоваться аналогией между влиянием дисперсии и реальной вязкости. Так, уменьшение амплитуды плоской бегущей волны из-за процессов вязкой диссипации определяется соотношением [8]:
где а — коэффициент поглощения. Величину этого коэффициента нетрудно рассчитать для конкретных значений параметра S и расстояния х1 по данным на рис. 2 и затем использовать для оценки величины Др на других расстояниях. Так, по материалам расчетов для координаты Х1 = 30 была получена следующая зависимость:
Примеры прогноза амплитуды сигнала на расстоянии Х1 = 60 с использованием (2.2) для значений ^ S = -2.3, -2, -1.865 и -1.5 показаны крестиками, которые практически точно легли на сплошные кривые.
После определения условий существенного влияния дисперсионных эффектов естественно рассмотреть вопрос о влиянии вязкости и теплопроводности реального газа на амплитуду бегущей волны. Соответствующие расчеты для Х1 = 60 представлены кривой 3 в зависимости от \% Яе в масштабе левой оси ординат. Расчеты выполнены при Рг = 0.7, Др1 = 10-3 и значении ^ S = -2.9, которое в приближении уравнений Эйлера обеспечивает воспроизведение амплитуды звуковой волны в конце канала точностью ~ 0.2%. При этом крестиками, как и в предыдущем случае, нанесены результаты прогноза на основе соотношения (2.1) величины вязкого затухания в конце канала для значений ^ Яе = 4.25, 4.5 и 4.75. Здесь для расчета коэффициента поглощения была использована формула Стокса — Кирхгоффа (2.13) из [8], в которую согласно рассматриваемой модели течения добавлен сомножитель X, что позволяет представить ее в следующем виде:
Др = Др1ехр(-ах1),
(2.1)
(2.2)
).
(2.3)
2.2
1.4
1.8
-2.6
-2.2
-1.8
1йДр
Рис. 4
Поскольку соотношения (2.2) — (2.3) для коэффициентов поглощения не зависят от амплитуды сигнала, они пригодны для оценки влияния как дисперсионных эффектов, так и физической вязкости при распространении сигналов с другой начальной амплитудой, но не превышающей значительно приведенную выше. Подобное ограничение обусловлено следующим обстоятельством.
Характерное изменение формы сигнала в выходном сечении канала по сравнению с исходным, о котором свидетельствует кривая 4 на рис. 1, связано, как известно, с более высокими температурой газа и скоростью звука в гребнях звуковой волны по отношению к минимумам. При дальнейшем распространении волны подобное нарастание крутизны ее переднего фронта должно привести к возникновению слабого разрыва. Формирование такого разрыва в виде локального пика давления на переднем фронте волны показано в верхней части рис. 3 кривой 1. Расчеты выполнены при ^ = -3, Др1 = 5.8 • 10- для невязкого газа. Ниже
показаны характерные этапы последующей эволюции формы волны, которые качественно до деталей согласуются с представленными в [8] на рис. 3.5 фотографиями осциллограмм сигнала, полученных экспериментально при возбуждении плоской ультразвуковой волны в воде. Такое соответствие свидетельствует об общности нелинейных явлений, определяющих изменение формы сигнала при его распространении независимо от характеристик конкретной среды.
В качестве критерия перехода звуковой волны в разрывное состояние можно использовать момент, когда амплитуда локального пика давления на ее переднем фронте начинает превосходить амплитуду основного сигнала. В частности, для кривой 1 этот момент соответствует координате х* = 32.5. Используя подобный критерий, можно провести
параметрические расчеты и определить характерные расстояния, на которых происходит качественное изменение формы сигнала, как функцию его исходной амплитуды. Результаты соответствующих расчетов в логарифмических переменных ^ Др1, ^ х* показаны на рис. 4 прямыми крестиками. Здесь же штриховой линией нанесена прямая
которая следует из аналитического решения (1.19), полученного в [8] для расстояния образования разрыва, если его представить в безразмерном виде. Хорошее соответствие между численными данными и аналитическим решением (2.4) свидетельствует о пригодности сформулированного выше критерия для оценки результатов численного моделирования.
(2.4)
1
0
1
51
О
-1
1
г2
О -1
О 2 4 6 К Х2
Рис. 5
Практическое использование как соотношения (2.4), так и его графического представления на рис. 4, зависит от степени влияния вязкой диссипации. Так, в верхней части рис. 3 кривыми 2 и 3 со сдвигом по координате Х1 в 0.2 и 0.4 представлены результаты расчета формирования разрыва в вязком теплопроводном газе при ^ S = -3, Ар1 = 5.8 • 10-3, Pr = 0.7 и числах Re = 106 и
3 • 105 соответственно. В первом случае влияние диссипации на амплитуду локального пика давления приводит к смещению точки формирования разрыва на расстояние в длину волны. Во втором случае подобное смещение составляет две длины волны. При Re = 107 воспроизводится состояние невязкого газа. Поскольку представленные материалы соответствуют сигналам с нулевой начальной фазой, то погрешность определения точки формирования разрыва соответствует длине волны. В связи с этим влиянием вязкой диссипации на формирование разрыва можно пренебречь, если это влияние не превышает длину волны. В свою очередь, число Re, при котором влияние диссипации соответствует длине волны, естественно выбрать в качестве граничного. Параметрические расчеты в диапазоне Др1 = (0.17 1.7) • 10-2 позволили определить
граничные значения чисел Яе*, которые в координатах ^ Др1, ^Яе* нанесены на рис. 4 косыми крестиками. Эти точки хорошо ложатся на прямую
^Яе* = 1.81 — 1.88 ^ Др1,
показанную сплошной линией.
Представленные выше материалы в определенной степени пригодны для оценки состояния звуковой волны при ее распространении по каналу переменного сечения. В качестве примера рассмотрим канал типа сопла Лаваля, основные пропорции которого соответствуют конфигурации плоского сопла трансзвуковой аэродинамической трубы Т-128 ЦАГИ, формирующего рабочий поток с числом М = 1.47 [5, 6]. Полуширина этого канала показана в верхней части рис. 5 как функция координаты Х2, нормированной на величину критического сечения 5*, которая в данном случае является характерным линейным масштабом. При этом входное сечение имеет размер 3 5*, выходное — 1.17 5*, положение критического сечения
соответствует Х2 = 5.5, а форма контура аппроксимирована аналитическими функциями. Ниже сплошной линией 1 представлены результаты расчета распространения звуковой волны в неподвижном невязком газе, заполняющем такой канал, при параметрах Др1 = 10-2, / = 1.43 и ^ S = -2.444. Качество этого расчета можно оценить на основе закона сохранения звуковой энергии в лучевой трубке (3.43) из [3], который
в одномерном приближении будет иметь вид:
(1 + -)2 = const. (2.5)
2рс с
Согласно (2.5) при и = 0 и рс = const амплитуда сигнала в критическом сечении должна возрасти в 1.73 раза. Линия 1 показывает, что именно это значение и воспроизводится. Далее, в расширяющейся части канала, при Х2 = 12.5X на переднем фронте волны формируется слабый разрыв, а амплитуда ее отрицательной части снижается до значения 1.59Ap1, которое на 0.6% меньше прогнозируемого 1.60Api. Данное расхождение следует отнести как на счет диссипации энергии волны в разрыве, так и на счет дисперсионных эффектов. Так, для канала постоянного сечения при lg S = -2.444 и х2 = 12X по соотношениям (2.2) и (2.1) получим а1 = 2.6•lO-4 и Si = 0.997, т. е. влияние дисперсии составляет половину наблюдаемого уменьшения амплитуды. В свою очередь, значению амплитуды Ap1 = 1.59 • 10-2 по соотношению (2.4) соответствует расстояние формирования разрыва х* = 11.7X, которое в пределах существующей погрешности совпадает с указанным выше.
В заключение рассмотрим пример численного моделирования распространения звуковой волны от сечения при х = 0 в дозвуковой части вниз по потоку, который реализуется в данном сопле при соответствующем перепаде давления. Расчет проводился в два этапа. На первом этапе методом установления рассчитывался стационарный трансзвуковой поток с использованием граничных условий (1.3). Установление доводилось до состояния постоянства расхода и сохранения полной энергии с точностью 10-6. Затем в нулевом узле граничные условия формировались согласно (1.4).
Так, в нижней части рис. 5 сплошной линией 2 представлены результаты расчета распространения звуковой волны в ускоряющемся трансзвуковом потоке при Ap1 = 10-2, f = 1.43 и lg S = -2.52. Значения сигнала в конкретной точке S2 представляют собой нормированную на Ap1 разность между мгновенным локальным значением безразмерного статического давления p(x2, t) и его значением в данной точке при стационарном течении p(x2, 0), распределение которого показано сплошной линией 3. При этом в сверхзвуковой части канала длина волны относительно
неподвижной системы координат возрастает практически вдвое за счет сложения локальной скорости звука и такой же скорости потока. В результате в пределах сопла реализуется девять периодов волны и форма сигнала сохраняется близкой к гармонической. Это позволяет воспользоваться соотношением (2.5) и оценить степень сохранения потока средней звуковой энергии
в процессе распространения сигнала вниз по течению. Расчеты, выполненные в диапазоне f = 0.36 ^ 1.44, Ap1 = (0.1 ^ 1.73) • 10-2 при lg S = -3.52 ^ -2.92, показали, что поток средней звуковой энергии в критическом и выходном сечениях по отношению к потоку во входном сечении сохраняется с относительной погрешностью ~ 10-2, которая согласуется с общей точностью проведения расчетов.
Если сверхзвуковую часть сопла продолжить каналом постоянного сечения достаточной протяженности, можно достичь состояния формирования разрыва и при наличии потока, о чем свидетельствует кривая 4. В свою очередь, при соответствующем увеличении частоты переход к разрывной форме сигнала будет осуществляться в пределах сопла. Более детальный анализ подобных состояний предполагается осуществить при численном исследовании взаимодействия звуковых колебаний и спонтанной конденсации трансзвукового потока в сопле трубы Т-128 ЦАГИ.
Работа выполнена при поддержке РФФИ (проекты № 04-01-00810-а и № 05-08-33663-а).
ЛИТЕРАТУРА
1. Корценштейн Н. М. Анализ результатов численного моделирования конденсационной релаксации пересыщенного пара // Коллоидный журнал. — 2002. Т. 64,
№ 5.
2. Блохинцев Д. И. Акустика неоднородной движущейся среды. — 2-е изд. —
М.: Наука. —1981.
3. Осташев В. Е. Распространение звука в движущихся средах. — М.: Наука. — 1992.
4. Завадский В. Ю. Метод конечных разностей в волновых задачах акустики. — М.: Наука. — 19S2.
5. Verhovskij V. P., Filipenkov V. N., Chirihin A. V. Nonequilibrium steam condensation in the humid air flow in the flat nozzle of high dimensional transonic tunnel. Fundamental Research in Aerospace Science // International Conference. Book Abstracts. Section 3: Nonequilibrium Flows and Rarefied Gas Flows. TsAGI. — 1994.
6. Чирихин А. В. Неравновесные течения влажного запыленного воздуха в сопле крупномасштабной трансзвуковой аэродинамической трубы// Изв. РАН, МЖГ. — 1999, № 4.
7. Крайко А. Н., Осипов А. А. Исследование отражения возмущений от дозвуковой части сопла Лаваля // Изв. РАН, МЖГ. — 1973, № 1.
S. Красильников В. А., Крылов В. В. Введение в физическую акустику. — М.: Наука. — 19S4.
9. Лойцянский Л. Г. Механика жидкости и газа. — М.: Наука. —1970.
10. Рычков А. Д. Математическое моделирование газодинамических процессов в каналах и соплах. — Новосибирск: Наука. —19SS.
11. Березин И. С., Жидков Н. П. Методы вычислений. — 2-е изд. — М.: Гос. изд. физ.-мат. лит. — 1962. Т. 1.
12. Лаваль П. Нестационарный метод расчета трансзвуковых течений в соплах. Численные методы в механике жидкостей. — М.: Мир. — 1973.
Рукопись поступила J8/XI2005 г.
