CC BY
58
Вычислительные технологии
Том 17, № 2, 2012
Численное моделирование течения газа
>к
около пластины с микровдувом*
А. В. Блзовкин, В.М. Ковеня, А. С. Лебедев Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: eli-ot@yandex.ru, kovenya@ict.nsc.ru, sasa@ict.nsc.ru
Численно моделируются течения газа около пластины при малых скоростях. Исследуется влияние интенсивности микровдува с части поверхности и его распределение на поведение локального и интегрального коэффициентов трения. Расчеты показали, что при сохранении общей интенсивности вдува интегральное значение трения практически не зависит от распределения участков вдува по поверхности. Полученное в расчетах снижение сопротивления трения качественно и количественно согласуется с результатами эксперимента.
Ключевые слова: уравнения Навье — Стокса, численные методы, моделирование течений газа, микровдув, сопротивление трения.
Введение
Поиск путей уменьшения аэродинамического сопротивления обтекаемых тел остается одной из актуальных проблем современной прикладной гидро- и аэродинамики. К настоящему времени разработан и апробирован целый ряд эффективных методов управления течением для уменьшения сопротивления трения. В последние десятилетия определенный прогресс в этом направлении был достигнут посредством подходов, включающих различную технику управления турбулентными пограничными слоями: использование механизма периодического вдува/отсоса, устройств разрушения вихрей, макропористых стенок и других пассивных и активных методов с целью модификации пограничного слоя для снижения сопротивления трения обтекаемых тел [1-12]. В работах [5-7] предложен метод управления пограничным слоем путем вдува через микропористую проницаемую поверхность, имеющую высококачественную обработку. Экспериментальные исследования, выполненные в [5-10], показали, что поверхностное трение можно снизить почти на 50 %, и тем самым продемонстрировали высокую эффективность данного подхода. Отметим, что в приближении полных уравнений Навье — Стокса данный класс задач до настоящего времени моделировался лишь в рамках модели сжимаемого газа [8, 12].
В настоящей работе проведено численное моделирование течения газа около пластины при малых скоростях. Цель моделирования — исследование влияния интенсивности микровдува с части поверхности и характера его распределения на поведение локального и интегрального коэффициентов трения. Известно, что в классической постановке такие течения изучаются в приближении пограничного слоя. Однако в зависимости от
* Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00294-а) и Интеграционного проекта СО РАН № 130.
интенсивности вдува возможны режимы течения, в которых приближение пограничного слоя может стать несправедливым, поэтому в данной работе в качестве математической модели выбраны уравнения Навье — Стокса вязкой несжимаемой жидкости. Сравнение результатов численных расчетов по моделям Навье — Стокса и пограничного слоя показало, что при увеличении расхода микровдува необходимо использовать полные уравнения Навье — Стокса вязкой несжимаемой жидкости. Проведенные расчеты течений в рамках модели Навье — Стокса вязкой несжимаемой жидкости около тел с микровдувом газа с части поверхности показали качественное и количественное совпадение с результатами экспериментов.
1. Постановка задачи и метод решения
Задача об обтекании пластины потоком газа с микровдувом с части поверхности рассматривалась в рамках уравнений пограничного слоя и полных уравнений Навье — Стокса вязкой несжимаемой жидкости, которые в двумерном случае записываются следующим образом:
ЗУ 1 ЗУ 2
дХГ + = 0-
Зу1 д , 2 , д
+ дХ1 (у2 + р - Т11) + дХ2(У1У2 - Т12) = 0'
дУ2 д д , 2 ,
+ дХ1(У1У2 - Т21) + + р - Т2^ =0' (1)
(дУ дуД
где у — компоненты скорости; р — давление; тц = я —--+ —— — тензор напряже-
\ дх] дХг I
ний; я = (я + Ят)/Ке — коэффициент вязкости (сумма ламинарной и турбулентной вязкости); Ие= — число Рейнольдса (Ц — скорость набегающего потока, Ь —
характерный размер области). В качестве характерного размера Ь выбран 1 м. Плотность полагалась постоянной и из уравнений исключена. При переходе к безразмерному виду все величины нормировались на их невозмущенные значения.
В двумерном приближении расчетная область имела форму прямоугольника, основанием которого служила пластина конечной длины. На части поверхности этой пластины задавался микровдув газа (рис. 1). Параметры течения во входном сечении брались из решения двумерных уравнений пограничного слоя [13]. На верхней границе области, расположенной вне пограничного слоя, задавались условия невозмущённого потока. В выходном сечении для вектора скорости задавались "мягкие" условия [14]. На поверхности пластины, кроме части области, из которой осуществлялся вдув, задавались условия прилипания. На участке вдува использовалось граничное условие
У1 =0, У2 = Сь.
Интенсивность вдува Сь представляет собой отношение скорости вдуваемого газа к скорости набегающего потока. В расчетах величина Сь не зависела от времени и изменялась в пределах 0 < Сь < 10-2, что включало диапазон значений вдува в экспериментах [8].
Положение верхней границы менялось в зависимости от интенсивности микровду-ва. Тестовые расчеты показали, что для исключения влияния верхней границы на решение, получаемое вблизи пластины, необходимо располагать эту границу на расстоя-
Рис. 2. Поведение местного коэффициента трения при различном положении верхней границы расчетной области: кривые 1, 2, 3, 4 относятся к случаям, когда граница находилась на расстоянии 2.5, 5, 10 и 20 толщин пограничного слоя соответственно
числялся по формуле Cf = 2^
нии 10-20 толщин пограничного слоя. В качестве примера на рис. 2 показано поведение местного коэффициента трения на участке пластины в области вдува при интенсивности вдува Сь = 1.5 х 10-3 (рис. 2, а) и Сь = 3.0 х 10-3 (рис. 2, б). Течение в данном примере ламинарное, число Рейнольдса Яе = 105. Местный коэффициент трения вы-
ЗУЛ
— . Кривая 1 соответствует расчету, в котором
х2
/ Х2 =0
верхняя граница расположена на расстоянии около 2.5 толщин пограничного слоя, с возрастанием номера кривой это расстояние удваивается, т. е. для кривой 2 оно равно 5 толщин пограничного слоя и т. д. Из рисунка видно, что кривые 3 и 4 близки, поэтому отодвигать верхнюю границу на ещё большие расстояния не имеет практического смысла.
Численное решение уравнений (1) находилось конечно-разностными методами. Наибольшие градиенты решения имеют место в тонком слое вблизи пластины, поэтому для
получения достаточной точности сетка возле пластины сгущалась. С этой целью вводилось преобразование координат, переводящее расчетную область в единичный квадрат. Уравнения решались в преобразованных координатах. В вертикальном направлении преобразование координат (сгущение) задавалось формулой
Х2 (k) = X2min + (X2 max - X2 min )(qfc/fcmax - 1)/(q - 1), q > 1,
где (X2max — X2min) — высота расчётной области в исходных координатах. Данная формула обеспечивает сгущение сетки по нормали к пластине, причем шаг сетки у поверхности пластины приблизительно в q раз мельче шага у верхней границы. Параметры сгущения выбирались таким образом, чтобы в область пограничного слоя попадало около половины узлов. В продольном направлении x1 сетка была либо равномерной, либо со сгущением в области вдува.
Для обеспечения большей достоверности расчётов решение уравнений (1) находилось по двум численным алгоритмам. Оба алгоритма аппроксимируют исходные уравнения со вторым порядком по пространственным координатам. Первый алгоритм основан на методе приближенной факторизации при специальном расщеплении исходных уравнений по физическим процессам и пространственным направлениям (подробнее см. [14]). Во втором алгоритме, использующем разнесенную сетку и идею искусственной сжимаемости, давление на новом временном слое и компоненты скорости, удовлетворяющие уравнению неразрывности, отыскиваются с помощью внутренних итераций.
Суть второго алгоритма поясним, записав двумерные уравнения Навье — Стокса следующим образом:
dvi dv2
&Г + eü = "•
öv1 öp
öt öx1 ь
— + = R2. (2)
öt öx2
В уравнениях (2) искомыми величинами являются давление p и компоненты скорости v1 и v2 вдоль пластины и по нормали к пластине. В уравнениях движения (второе и третье уравнения) вязкие и конвективные слагаемые находятся в правой части и обозначены R1 и R2. Считаем, что в расчетной области введена сетка, узлы которой имеют координаты x1i, x2*. Давление определим в центре ячейки, а скорости — на соответствующих гранях, т.е. искомыми будут сеточные функции pi+1/2,fc+1/2, (^1)^+1/2, (^2)^+1/2,*. Если эти функции известны в момент tn = тп, где т — шаг дискретизации по времени, то в момент tn+1 = tn + т они находятся из системы уравнений
—1--1--2— = 0,
öx1 öx2
vn+1 - < öpn+1 Rn т öx1 1'
vn+1 - vn + =Rn. (3)
т öx2
Уравнение неразрывности в (3) аппроксимируется в центре ячейки, а уравнения движения — на соответствующих гранях. Все пространственные производные аппроксимируются со вторым порядком. Величины рП+1/2 &+1/2, (У1)Г'++1/2, )П+1/2 к, удовлетворяющие (3), находятся с помощью итераций из уравнений
р^1 - рГ , д<+1 д^1 п
ат дх1 дх2 '
уГ+1 - у? др/
1_1 | ^ = дп
т дх1 1'
- у? др/
V2-2 + = (4)
т дх2
По известному р^ (в качестве начального значения берется р0 = рп) из уравнений движения в (4) находятся и У2+1, затем из уравнения неразрывности находится давление р^+1 и т. д. Итерации по индексу V выполняются до тех пор, пока с заданной точностью не будет удовлетворяться уравнение неразрывности.
Для проверки точности использованных алгоритмов в диапазоне чисел Рейнольдса 103 < И,е < 106 были проведены расчеты обтекания пластины на сетках, обеспечивающих выполнение закона сохранения массы с точностью до 10-4 -10-6 (как правило, было достаточно задать 100 х 100 узлов на неравномерной сетке).
2. Результаты численных расчетов 2.1. Ламинарный режим
Расчёты обтекания пластины для ламинарного режима течения при отсутствии вдува, выполненные на основе двумерных моделей пограничного слоя и Навье — Стокса, показали практически полное совпадение результатов во всем диапазоне чисел И,е. При наличии вдува появляются различия, возрастающие с увеличением интенсивности вдува. Для иллюстрации на рис. 3 приведено распределение местного коэффициента трения
0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Х^ 0.4 0.6 0.8 1 1.2 1.4 1,6 1.8 Л
Рис. 3. Распределение коэффициента трения на пластине при различных значениях Сь: 1 — без вдува, 2 — модель Навье — Стокса, 3 — модель пограничного слоя
для Сь = 1.0 х 10-3 и Сь = 3.0 х 10-3 при И,е = 105. Область вдува находилась на участке 0.7 < х1 < 0.91 и на рисунке обозначена вертикальными линиями. В модели Навье — Стокса в отличие от модели пограничного слоя влияние вдува распространяется вверх по потоку, что приводит и к различию значений коэффициента трения. При небольших интенсивностях вдува интегральный коэффициент трения, вычисленный по модели пограничного слоя, оказывается близок к значению, полученному из решения уравнений Навье — Стокса. Например, при Сь = 1.0 х 10-3 значения коэффициента трения по мере удаления от участка вдува различаются всё меньше, а в окрестности области вдува (на участке 0.4 < Х1 < 1.4) отличие интегрального коэффициента трения составляет менее 0.5 %. Возрастание интенсивности вдува приводит к увеличению различий в поведении коэффициента трения и несовпадению интегрального коэффициента трения, что объясняется неприменимостью модели пограничного слоя для течений со вдувом газа с поверхности. В расчётах с другим размером и положением области вдува (вдув осуществлялся с 55.6 % части поверхности) были получены аналогичные результаты.
Для сравнения коэффициентов трения при различной организации вдува были проведены расчеты, в которых в первом случае вдув задавался равномерным по всей подобласти, во втором — осуществлялся из последовательно расположенных в области вдува участков, количество которых варьировалось, но общая площадь не менялась и составляла 19 % от площади равномерно распределенного вдува. При изменении пространственной структуры изменялась и интенсивность вдуваемого газа таким образом, чтобы массовый расход сохранился. Область вдува (как равномерного со всего участка, так и из отдельных участков) располагалась на участке 0.7 < Х1 < 0.91.
На рис. 4 показано распределение коэффициента трения с/ для ламинарного режима течения при И,е= 105 и Сь = 10-3. Кривая 1 соответствует течению без вдува, 2 — равномерному вдуву, кривые 3 и 4 относятся к областям вдува из 4 и 16 участков. При расчете течения с несколькими подобластями вдува использовались более мелкие сетки при условии, чтобы на каждую подобласть вдува приходилось от 4 до 8 узлов (ячеек). Этого оказалось достаточно, чтобы в зоне вдува распределение коэффициента трения не зависело от числа узлов сетки в зоне вдува. Расчёты показали, что с увеличением числа участков вдува (при сохранении массового расхода) интегральный коэффициент трения приближается к соответствующему значению при равномерном вдуве.
Ниже приведено относительное отклонение 8 = (С/ — Cfu)/Cfu интегрального коэффициента трения пластины С/ при разном количестве участков вдува от коэффициента С/и для равномерного вдува при И,е= 105 и Сь = 10-3 (ламинарный режим). Местный
2.8 г
ш
0.6
0.7
0.8
Й.,.9
1.0
1.1
Рис. 4. Распределение коэффициента трения с/ при разном числе участков вдува
коэффициент трения в расчетах различался на интервале 0.6 < х1 < 1.2, поэтому инте-
1.2
гральный коэффициент трения вычислялся как С/ = J С/
0.6
Число участков вдува 4 16
5 = (С/ - С/„)/С/„ -0.45 % -0.35 %
Из этих данных видно, что если область вдува состоит из большого числа участков, то коэффициент сопротивления приближается к значению, соответствующему равномерному вдуву при условии сохранения массового расхода вдуваемого газа. Полученные результаты позволяют сделать вывод о правомочности замены дискретно распределенного вдува (при большом числе участков) на равномерный вдув в данной подобласти. Это упрощает решение задачи, так как дает возможность использовать меньшее число узлов сетки при сохранении точности расчетов.
2.2. Турбулентный режим
В расчетах, направленных на исследование влияния турбулентности на течение газа около пластины со вдувом с части поверхности, использовалась алгебраическая модель турбулентности — модель пути перемешивания Прандтля [13, 15]. В исходных уравнениях к вязкости добавляется турбулентная вязкость вида = /2|дг'1/дж21, где I — путь перемешивания, во внешней области равный
те
I = 1.1875 У ^(1 -
0
а в пристеночной области —
/ = 0.4x2 (1 - е-в) , в =26
\
Зг>1
Ие 0.75
+ 0.25— дх
Х2=0
Переход от одной формулы к другой осуществлялся из условия непрерывности турбулентной вязкости. На рис. 5 показано распределение коэффициента трения при турбулентном обтекании пластины для разного числа участков вдува при Ие= 105 и Сь = 10-3 (здесь линия 1 — течение без вдува, 2 — равномерный вдув, линии 3, 4 — дискретный вдув из четырех и восьми участков).
Ниже приведено относительное отклонение 8 = (С/ - С/и)/С/и интегрального коэффициента трения для четырех, восьми и шестнадцати участков вдува (турбулентный режим):
Число участков вдува 4 8 16
5 = (С/ - С/и)/С/„ 0.09 % -0.04 % -0.04 %
Из этих данных, а также из рис. 5 следует: в случае турбулентного течения интегральный коэффициент поверхностного трения с ростом числа участков вдува стремится к его значению при равномерном вдуве, что позволяет в дальнейшем ограничиться расчетами с равномерно распределенным вдувом.
В следующих расчетах исследовано влияние интенсивности вдува на коэффициент поверхностного трения. Вдув осуществлялся с участка 0.7 < х1 < 0.91. На рис. 6 приведен местный коэффициент трения с/ для различных интенсивностей вдува Сь: 1 — 0, 2 — 1.0х10-3, 3 — 1.5х10-3, 4 — 3.0х10-3, 5 — 5.1 х10-3, 6 — 1.0х10-2 при Яе= 105. Кривые I соответствуют расчетам по методу искусственной сжимаемости, II — по методу расщепления. Наблюдается хорошее совпадение результатов, полученных различными алгоритмами.
На рис. 7 сравнивается местный коэффициент трения, полученный в расчетах с использованием уравнений Навье — Стокса (I) и пограничного слоя (II), для различных интенсивностей вдува Сь: 1 — 0, 2 — 1.0 х 10-3, 3 — 1.5 х 10-3, 4 — 3.0 х 10-3. Видно, что с ростом интенсивности вдува различие значений коэффициента трения, рассчитанного по разным моделям, возрастает, причем модель пограничного слоя дает меньшее
Рис. 5. Поведение коэффициента трения с/ при разном числе участков вдува
/ — 1
\ б / —
1
к-.:"
Г ---
— I
— II
0.4
0.6
0.8
1.0
1.2
1.4
1.6
0.4
0.6
0.8
1.0
1.2
1.4
1.6
Рис. 6. Коэффициент поверхностного трения для различных интенсивностей вдува, рассчитанный по разным численным алгоритмам
Рис. 7. Коэффициент поверхностного трения для различных интенсивностей вдува, рассчитанный по разным математическим моделям
снижение коэффициента трения (для ламинарных течений модель пограничного слоя предсказывала большее снижение).
В следующей серии расчётов было проведено сравнение численных решений уравнений Навье — Стокса с результатами натурного эксперимента [8]. Начало пористой вставки, из которой осуществлялся поверхностный вдув, находилось на расстоянии 1181.3 мм от передней кромки пластины. Длина пористой вставки составляла 400 мм, а общая длина пластины — 2313.6 мм. Эксперименты проводились при скорости невозмущенного потока U= 21.2 м/c, что соответствовало числу Рейнольдса 1.4 х 106. При отсутствии вдува коэффициенты трения, полученные на основе решений уравнений Навье — Стокса и уравнений пограничного слоя, совпали. На рис. 8 показан профиль скорости при Rex = 2.1 х 106, полученный из решения уравнений Навье — Стокса (кривая 1), представленный в универсальных переменных закона стенки u+ = f (log y+), где u+ и y+ в безразмерных переменных имеют вид
u
+
V\
\
1 dvi
Re dx2
У+ =
X2
X2=0
\
Re
dvi
dx2
X2=0
Для сравнения на рис. 8 также приведены профиль и+ = у+, соответствующий линейному распределению скорости в вязком подслое (кривая 2), и профиль и+ = 5.75у+ + 5.5 [15], отвечающий логарифмическому участку пограничного слоя (кривая 3). Видно, что профиль скорости, полученный в расчетах по модели Навье — Стокса, находится в хорошем соответствии с классическими данными, что иллюстрирует достаточно высокую точность численных результатов.
На рис. 9 представлено изменение коэффициента трения Cf для режима течения без вдува (1-3) и со вдувом (4-7), интенсивность которого равна Сь = 1.48 х 10-3, что в эксперименте соответствовало массовому расходу Q = 0.038 кг/(с • м2). В размерных физических величинах удельный массовый расход определяется как Q = ри^у2. Продольная координата нормирована следующим образом: X = (Х\ — 0.1091)/— 0.1091),
Рис. 8. Профиль скорости в переменных закона стенки
Рис. 9. Изменение коэффициента трения Cf по Рис. 10. Зависимость коэффициента трения Cf длине модели от расхода (
где Ь = 2.3136 — длина пластины (область вдува выделена вертикальными пунктирными линиями). Здесь данные эксперимента — кривые 1, 4 [8], расчёты по модели Навье — Стокса методом искусственной сжимаемости — кривые 2, 5, методом расщепления — кривые 3, 6 и расчёты по модели пограничного слоя — кривые 3, 7. Отличие рассчитанного коэффициента трения Cf от экспериментального составляет около 7%.
На рис. 10 показано изменение коэффициента трения в середине пористой вставки в зависимости от интенсивности поверхностного вдува (здесь 1 — расчёт по модели пограничного слоя, 2 — по модели Навье — Стокса, 3 — эксперимент). Отметим, что с ростом интенсивности вдува результаты расчетов по модели Навье — Стокса приближаются к экспериментальным данным.
Заключение
Численное моделирование турбулентных течений со вдувом газа с части поверхности в приближении уравнений Навье — Стокса вязкой несжимаемой жидкости даёт хорошо согласующиеся с экспериментом результаты, что позволяет утверждать, что снижение сопротивления трения обтекаемых поверхностей с помощью микровдува является эффективным способом управления турбулентным пограничным слоем.
Проведенные расчеты и сравнения с экспериментом показали, что применимость модели пограничного слоя ограничена относительно небольшими интенсивностями вду-ва. Для ламинарных режимов течения с увеличением интенсивности вдува эта модель предсказывает большее снижение коэффициента трения, чем модель Навье — Стокса, а для турбулентных — наоборот, дает заниженное по сравнению с моделью Навье — Стокса (и с экспериментом) снижение коэффициента.
Влияние микровдува на коэффициент трения пластины на участке ниже по течению от области, с которой осуществляется микровдув, является более сильным, чем на участке, расположенном выше по потоку, причем влияние возрастает при увеличении интенсивности вдува. Данный факт может быть использован для создания более эффективного метода снижения турбулентного трения.
Авторы выражают благодарность д.ф.-м.н. В.И. Корнилову за предоставленные экспериментальные данные и к.ф.-м.н. А.Н. Попкову за результаты расчетов в рамках
модели пограничного слоя.
Список литературы
[1] Корнилов В.И. Проблемы снижения турбулентного трения активными и пассивными методами (обзор) // Теплофизика и аэромеханика. 2005. Т. 12, № 2. C. 183-208.
[2] Kornilov V.I., Boiko A.V. Periodic forcing of the turbulent boundary layer on a body of revolution // AIAA J. 2008. Vol. 46, Ко. 3. P. 653-663.
[3] Бойко А.В., Корнилов В.И. Влияние периодического вдува/отсоса через последовательно расположенные кольцевые щели на турбулентный пограничный слой тела вращения // Теплофизика и аэромеханика. 2008. Т. 15, № 1. C. 11-29.
[4] Шиплюк А.Н., Буров Е.В., Млслов А.А., Фомин В.М. Влияние пористых покрытий на устойчивость гиперзвуковых пограничных слоев // ПМТФ. 2004. Т. 45, № 2. С. 169-176.
[5] Hwang D.P. A Proof of Concept Experiment for Reducing Skin Friction by Using a Micro Blowing Technique. AIAA Paper. 97-0546. 1997.
[6] Hwang D.P. An Experimental Study of Turbulent Skin Friction Reduction in Supersonic Flow Using a Microblowing Technique. AIAA Paper. 2000-0545. 2000.
[7] Hwang D.P. Review of Research into the Concept of Micro-Blowing Technique for Turbulent Skin Friction Reduction. Invited Paper, Prog. in Aerospace Sciences. October, 2004.
[8] Корнилов В.И., Бойко А.В. Использование микровдува воздуха через пористую стенку для снижения трения на плоской пластине // Вестник НГУ. Физика. 2010. Т. 5, вып. 3. С. 38-44.
[9] Корнилов В.И., Бойко А.В., Попков А.Н. Влияние микровдува воздуха через проницаемую стенку на турбулентный пограничный слой // Там же. 2011. Т. 6, вып. 1. С. 77-83.
[10] Lin Y.L., Chyu M.K., Shih T.I.P. et al. Skin Friction Reduction Through Micro Blowing. AIAA Paper 98-0359. 1998.
[11] Kudriavtsev V., Braun M.J. Computational study of micro-blowing for shear force reduction // VP-Vol. 424-2. Computational Technologies for Fluid Thermal Structural/Chemical Systems with Industrial Applications. ASME. 2001. Vol. 2.
[12] Li Jian, Lee Chun-Hian, Jia Liping, Li Xuzhi. Numerical study on flow control by micro-blowing // Proc. of the 47th AIAA Aerospace Sciences Meeting. AIAA 2009-779. Orlando, Florida.
[13] Колобов Б.П., Попков А.Н. Моделирование двумерного пограничного слоя // Численные методы в механике жидкости и газа. Новосибирск: ИТПМ СО АН СССР, 1980.
[14] Блзовкин А.В., Вавилова О.М., Ковеня В.М. Метод факторизации для численного решения уравнений вязкой несжимаемой жидкости // Вычисл. технологии. 2009. Т. 4, № 2. С.13-31.
[15] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1978. 736 с.
Поступила в 'редакцию 23 января 2012 г., с доработки — 5 марта 2012 г.
