Явные разностные схемы с переменными шагами по времени в задаче вихреразрешающего моделирования течения несжимаемой жидкости над шероховатой поверхностью Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 11, часть 2, Специальный выпуск, 2006

ЯВНЫЕ РАЗНОСТНЫЕ СХЕМЫ С ПЕРЕМЕННЫМИ ШАГАМИ ПО ВРЕМЕНИ В ЗАДАЧЕ ВИХРЕРАЗРЕШАЮЩЕГО МОДЕЛИРОВАНИЯ ТЕЧЕНИЯ НЕСЖИМАЕМОЙ ЖИДКОСТИ НАД ШЕРОХОВАТОЙ ПОВЕРХНОСТЬЮ*

К. В. Ушаков

Институт, вычислительной математики РАН, Москва, Россия,

e-mail: nucrectOinm.ras.ru

A three-dimensional large-eddy simulation of an incompressible fluid flow over a rough surface has been conducted. The method of variable time steps with the most extensive stability domain along the imaginary axis (among the 1-4-step methods) is applied.

Движение нейтрально стратифицированной несжимаемой жидкости с хорошей точностью описывается системой уравнений Навье — Стокса. Однако в случае турбулентных течений эта система не поддается аналитическому решению. Численное решение оказывается очень чувствительным к пространственной и временной точности представления вихревых структур различных размеров. Поэтому прямое численное моделирование таких течений является очень дорогостоящим, объем требуемых вычислительных ресурсов оказывается пропорциональным третьей степени числа Рейнольдса [1]. Так, для прямого моделирования пограничного слоя атмосферы с типичным числом Re = 109 количества расчетных точек и временных слоев будут иметь порядки 1020 и 107 соответственно, что неприемлемо для современной вычислительной техники. Альтернативой может служить предложенный во второй половине XX в. метод вихреразрешающего моделирования. Он подразумевает явное описание крупных вихрей и параметризацию мелких с использованием осредненных по пространству значений скорости. Однако и после такого осреднения задача остается ресурсоемкой, что вызывает необходимость использования параллельных вычислительных систем. В связи с этим в задачах такого типа могут оказаться полезными явные разностные схемы, поскольку они в отличие от неявных обычно допускают эффективное распараллеливание, в том числе в расчетных областях сложной формы. Использование переменных шагов по времени дает возможность ослабить жесткие условия устойчивости, ограничивающие шаг в случае его постоянства.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 04-05-64898, № 05-05-64918, № 05-01-00582) и программы РАН “Теоретические проблемы современной математики” (проект “Оптимизация вычислительных алгоритмов решения задач математической физики”).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

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

а(х, £) = У С(х,х!)а(х!,£)б1:х!, (1)

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

дщ _______д__________д_ др 1 д2щ

дЬ дх^1^ дх^ дх„• Кедх^дх^’

£ = * <3>

Ту = щщ - щ щ. (4)

Первые два уравнения выражают закон изменения импульса частицы жидкости и условие

неразрывности, В соответствии с равенством (4) в правой части (2) нелинейное слагаемое

разделено на составляющую, зависящую только от фильтрованных компонент скорости, и составляющую, которая включает не описываемые явно в рамках данного подхода величины.

Согласно подтверждаемому опытом предположению Колмогорова [1], при достаточно больших числах Рейнольдса в спектре трехмерной изотропной турбулентности существует так называемый инерционный интервал, в котором кинетическая энергия передается (без производства или диссипации) от крупномасштабных движений к мелкомасштабным. Принятие этой гипотезы дает возможность записать модельное замыкание, связывающее тензор напряжений с фильтрованными компонентами скорости Щ. В основу настоящей работы положена модель [2], использующая смешанное замыкание

т“- = ( — 2СА.2\Б\8г^) + Щ — Пг щ')а. (5)

о и » 77 1 (9й~г дщ\

Здесь индекс а означает анизотропную часть тензора; = - —----Ь — тензор

2 \ дх2 дхг _

скоростей деформации; |5'| = у^2А = (Д^Д^Д^)1/3 — среднегеометрическое значение шага сетки. Изотропная часть тензора напряжений считается включенной

в величину фильтрованного давления р. Пространственное распределение коэффициента турбулентной вязкости СА2^! вычисляется в динамической процедуре, использующей дополнительную “тестовую” фильтралцю [3], Поскольку кинетическая энергия решения системы уравнений Навье — Стокса в пределе больших чисел Рейнольдса не меняется со временем, к используемой численной схеме мы будем предъявлять требование консервативности, Пространственная дискретизация уравнений осуществляется с помощью консервативной схемы четвертого порядка [4],

Применяемые в расчетах модели явные схемы используют переменные шаги по времени, Процедура выбора последовательности шагов основана на анализе области устойчивости метода, соответствующего взятой последовательности [5]. Пусть имеется задача Коши для системы линейных обыкновенных дифференциальных уравнений, полученная из исходной задачи математической физики с помощью линеаризации и дискретизации по пространству:

пи

_АИ + /, и\,=0 = и°, *<Е[0,Т]. (6)

аъ

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

пк+1 — ик

тк+1

—Аик, к = 0,1,...,^ - 1, (7)

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

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

пжх |Рм(г)! < 1, (8)

г£ор+Л

где

N

Рм(г)=П(1 - Тгг)’ (9)

г=1

а Бр+Л — множество собственных значений матрицы Л, которые обладают неотрицательными действительными частями. Легко видеть, что многочлен Рм(Л) представляет собой оператор перехода всего цикла шагов вида (7): им = Рм(Л)и0. Множество комплексных чисел г, таких что |Рм(г)| < 1, будем называть множеством устойчивости многочлена Рм, Таким образом, построение устойчивого вычислительного метода сводится к построению такого многочлена вида, (9), множество устойчивости которого содержит в себе часть Л

те [5] приводится способ построения класса многочленов Рм, позволяющих для определенных задач получить ускорение счета в О^) раз по сравнению с методом Эйлера с постоянным шагом. Их области устойчивости окружают отрезки прямых и гипербол. Особенностью систем уравнений Навье — Стокса и вихреразрешающего моделирования

Л

в качестве Рм выбран многочлен

Р4(г) = 0,4 (1г), (10)

где 04(Ъ) = 1 - Ъ + Ъ2/2 - Ъ3/6 + Ъ4/24,

В соответствии с (9) временные шаги метода являются величинами, обратными к корням многочлена Р^{£). Многочлен (), (/) имеет наибольший отрезок устойчивости вдоль мнимой оси 2у/2г\ среди многочленов вида 1 — £+а£2+Ь£3+с£4 [5]. Множество устой-

чивости многочлена 04 изображено на рис, 1, для сравнения справа приведено множество

Рис. 1. Множества устойчивости Q4(t) и многочлена, соответствующего схеме Мацуно.

устойчивости многочлена, соответствующего схеме Мацуно, Изменение параметра I позволяет расширить или сжать множество устойчивости многочлена Р4(£) (при 0 < I < 1 и I > 1 соответственно), что приводит к обратно пропорциональному изменению всех временных шагов цикла,.

Корпи Р^^) могут быть комплексными. Поскольку коэффициенты многочлена действительны, МОЖНО считать, ЧТО если Т2І-1 комплексно, ТО Т2і = Т2І—1 ■ Заметим, что произведение пары соответствующих скобок в (9) можно представить в виде

Тогда пару комплексно-сопряженных шагов метода (7) можно реализовать в действительной арифметике по формулам, описывающим шаги по методу Эйлера и поправку:

Расчеты проводились па сетке размером 160 х 80 х 48 узлов с помощью параллельной ЭВМ, Модель реализована с периодическими горизонтальными граничными условиями, шероховатой поверхностью на нижней границе расчетной области и нулевым потоком импульса или с такой же шероховатой поверхностью (случай “плоского канала”) на верхней границе, В качестве начальных условий взято однородное течение под небольшим углом к длинной стороне расчетной области с наложением случайного шума. Поскольку вихреразрешающая модель, как было указано выше, дает на выходе фильтрованное поле скорости, обработаем это поле после завершения счета обращенным оператором фильтрации. Такая процедура помогает более точно воспроизвести в модели результаты экспериментальных исследований.

Модель генерирует нестационарное течение сложной структуры, В качестве примера на рис, 2 приведены характеристики поля горизонтальной скорости в плоскости, лежащей вблизи середины канала: стрелками показаны величина и направление, цветом на левом

(1 — т2і-і ^)(1 — т2і %) = (1 — Кіг)2 — ,

где

Ук+1/2 = ик — ^к+1Аик, Ук+1 = Ук+1/2 — ^к+1АУк+1/2, ик+2 = Ук+1 + 7к+1^к+1 (АУк+1/2 — Аик).

40

20

35

30

10

15

25

0,0

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

1-16,0 -14,6 -13,1 -11,7 -10,3 -8,82 -7,3В -5,94

А

10

20

30

40

50

60

70

80

04

Рис. 2. Величина и направление поля горизонтальной скорости вблизи середины канала: завихренность (слева), спектр (справа).

рисунке — вертикальная компонента завихренности, на правом — двумерный спектр поля горизонтальной скорости. (На левом рисунке изображена подобласть 80 х 40, на правом — по осям и шкале цвета отложены десятичные логарифмы величин. Для наглядности из скорости вычтено ее среднее значение.) Интенсивность высокочастотных пульсаций продольной компоненты скорости быстро спадает вблизи правого края спектра. Такой недостаток является общим для конечно-разностных моделей турбулентности. Он вызван тем, что пульсации, связанные с ошибками аппроксимации нелинейных членов баланса импульса, не отражают реальную физику процесса и не поддерживаются за счет нелинейных взаимодействий [2]. Поэтому для качественного воспроизведения спектра пульсаций необходимо достаточно высокое пространственное разрешение. При этом должны сохраняться и физические размеры расчетной области, иначе может возникнуть положительная обратная связь между значениями продольной компоненты скорости на ее левой и правой границах.

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

(п)

где осреднение производится по горизонтальным плоскостям; к = 0.4 — постоянная Кар-

а б

Рис. 3. Безразмерный градиент (а) и профиль скорости (б) по результатам расчетов. Для сравнения приведен точный логарифмический профиль — верхняя кривая на рис. (б). Номер горизонтального уровня отложен по вертикальной оси.

мана; и* — осредненное значение скорое ти трения и* (ж, у):

кл/иЦх, у, Аг/2) + иЦх, у, Аг/2)

и*(х,у) =

н'^

2о

— параметр шероховатости, принятый в данной работе равным 0.1м. Близость безразмерного градиента к единице выражает логарифмический характер зависимости осред-неппой скорости от высоты и более наглядна, чем прямое изучение профиля скорости.

На рис. 3 представлены графики безразмерного градиента и профиля продольной компоненты скорости для нижней половины канала, полученные по данным расчетов методом, соответствующим многочлену (10). Результаты говорят о хорошем качестве воспроизведения в модели логарифмического закона для профиля скорости, в целом не уступающем приведенному в [6]. Численные тесты также показали, что используемый алгоритм интегрирования по времени практически не вносит ошибок в баланс кинетической энергии решения. Для более детального выбора серии временных шагов представляется важным получение дополнительных сведений о спектре решаемой системы дифференциальных уравнений модели.

Список литературы

[1] КОЛМОГОРОВ А.Н. Локальная структура турбулентности в несжимаемой жидкости при очень больших числах Рейнольдса // Докл. АН СССР. 1941. Т. 30, № 4. С. 99-102.

[2] Глазунов А.В. Моделирование нейтрально стратифицированного турбулентного потока воздуха над горизонтальной шероховатой поверхностью // Изв. РАН. Сер. Физика атмосферы и океана. 2006. Т. 42. С. 307-325.

[3] PortE-Agel F., Meneveau С., Parlange М.В. A scale-dependent dynamic model for large-eddy simulation: application to a neutral atmospheric boundary layer // J. Fluid Meeh. 2000. Vol. 415. P. 261-284.

[4] Morinishi Y., Lund T.S., Vasilyev O.V., Moin P. Fully conservative higher order finite difference schemes for incompressible flow // J. Comp. Phys. 1998. Vol. 143. P. 90-124.

[5] Бахвалов H.C., Кобельков Г.М., Кузнецов Ю.А. и др. Численные методы решения задач математической физики // Совр. пробл. вычисл. математики и мат. моделирования. М.: Наука, 2005. Т. 1.

[6] Andren A., Brown A.R., Graf J. ет al. Large-eddy simulation of neutrally stratified boundary layer: a comparison of four computer codes // Quart. J. Royal Met. Soc. 1994. Vol. 120. P. 1457-1484.

Поступила в редакцию 9 ноября 2006 г.