ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 7, с. 1255-1265
УДК 519.634
ЧИСЛЕННАЯ СХЕМА ДЛЯ МОДЕЛИРОВАНИЯ ТУРБУЛЕНТНЫХ ТЕЧЕНИЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ С ИСПОЛЬЗОВАНИЕМ ВИХРЕРАЗРЕШАЮЩИХ ПОДХОДОВ1
© 2015 г. А. С. Козелков, В. В. Курулин
(607188Саров, Нижегородская обл., пр-т Мира, 37, ФГУП "РФЯЦ-ВНИИЭФ") e-mail: askozelkov@mail.ru; kurulin@mail.ru Поступила в редакцию 30.10.2014 г.
Исследуются свойства существующих численных схем дискретизации конвективных потоков и предложена альтернативная схема, формулировка которой выполнена на базе диаграмм нормализованной переменной. Анализируются диссипативные свойства предложенной схемы на примере моделирования турбулентных течений несжимаемой жидкости с использованием вихреразрешающих моделей. Показано, что при использовании предложенной схемы не возникает численных осцилляций решения, а ее диссипативные свойства превосходят свойства аналогичных схем с центрально-разностной аппроксимацией. Для предложенной схемы проведена калибровка модели крупных вихрей, как для свободных течений, так и для пристеночных течений, показана сеточная сходимость при решении задачи течения в плоском канале с использованием прямого численного моделирования турбулентности. Библ. 22. Фиг. 8. Табл. 1.
Ключевые слова: турбулентность, вычислительная гидродинамика, численная схема, дисси-пативность, вырождение изотропной турбулентности.
DOI: 10.7868/S0044466915070091
ВВЕДЕНИЕ
Среди вихреразрешающих подходов к моделированию турбулентности, за исключением метода прямого численного моделирования (см. [1]) (DNS), наиболее распространенными являются: метод крупных вихрей (см. [2]) (LES), метод отсоединенных вихрей (см. [3]) (DES), а также различные варианты зонных RANS-LES подходов (см. [4]). Причем моделирование с использованием первого из них, метода LES, является основой и для двух остальных.
Метод LES основан на решении отфильтрованных уравнений Навье—Стокса, причем в качестве размера фильтра, как правило, используется характерный размер расчетной ячейки (см. [2]). В таком подходе вихревые структуры, превышающие размеры фильтра, разрешаются численно, а все вихревые структуры с меньшими размерами моделируются с помощью той или иной под-сеточной модели турбулентности. Влияние малых вихрей на основное течение осуществляется с помощью добавления к молекулярной вязкости среды дополнительной подсеточной вязкости.
По сравнению с турбулентной вязкостью, определяемой при использовании RANS моделей, подсеточная вязкость, как правило, имеет на порядок меньшие значения, сопоставимые со значениями схемной вязкости, возникающей из-за решения дифференциальных уравнений на разностных сетках с конечным уровнем разрешения (см. [2]). В этом состоит главная проблема метода LES — величина подсеточной вязкости, которую необходимо привнести в уравнение движения напрямую зависит от схемной вязкости используемой численной схемы — чем больше схемная вязкость, тем меньше должен быть уровень дополнительной подсеточной вязкости (см. [2]). Стоит отметить, что далеко не каждая численная схема вообще способна обеспечить приемлемую величину схемной вязкости (или иначе — численной диссипации), не превышающую реальную диффузию, вызванной наличием вихрей малых масштабов.
1) Работа выполнена при финансовой поддержке РФФИ (код проекта 13-07-12079).
Практика применения вихреразрешающих подходов на основе метода LES большей частью сводится к расчетам на многоблочных структурированных сетках с применением схем высокого порядка (см. [5], [6]). Использование сеток такого класса вызывает определенные проблемы при моделировании течений в промышленных конструкциях сложной геометрической конфигурации, где построение блочно-структурированной сетки возможно, но это неоправданно дорого и крайне неэффективно (см. [7]). При расчетах индустриальных задач удобнее использовать неструктурированные расчетные сетки, состоящие из многогранников произвольной формы с выделением призматических пограничных слоев, а также LES-областей с изотропным, преимущественно шестигранным разбиением (см. [8]). Однако использование сеток такого класса приводит к ряду трудностей, одна из которых заключается в том, что на неструктурированной сетке заметно сужается круг доступных схем для дискретизации конвективных потоков.
Исследование свойств основных схем дискретизации конвективных потоков на неструктурированных сетках можно найти в [9]—[12]. В [11], [12] проводится анализ наиболее популярных схем, результаты которого показывают, что необходимый уровень диссипации могут обеспечить только схемы, основанные на центральных разностях. Использование центрально-разностных схем приводит к численным осцилляциям (см. [13]), которые, как правило, развиваются в пограничном слое и вблизи границ вход/выход. Уменьшения осцилляций и увеличения устойчивости центрально-разностной схемы можно добиться ее смешением с противопоточной (см. [14]). Данный способ широко применим на практике, однако, в [11] показано, что даже внесение 10% противопоточности увеличивает численную диссипацию схемы примерно в 2 раза.
Другой способ решения проблемы представлен в [9]. В ней предложена центрально-разностная схема с добавлением искусственной вязкости. Показано, что данная схема приводит к достаточно хорошим результатам для свободных течений при использовании гексагональных сеток, проверка ее диссипативности в случае пристеночных течений не приводится. В [10] описывается центрально-разностная схема с автоматическим смешением с противопоточной схемой, причем фактор смешения строится таким образом, что центрально-разностная схема включается в тех областях, где она наиболее устойчива и необходима для точного разрешения крупномасштабных вихревых структур. В пристеночной области, а также в области, свободной от отрывных вихрей, используется противопоточная схема первого порядка.
Хорошим способом монотонизации схем высокого порядка является применение ограничителей. В [13] описывается схема GAMMA, основанная на центрально-разностной схеме, ограниченной по критерию локальной ограниченности (Convection Boundedness Criterion, CBC (см. [15]) с применением диаграмм нормализованной переменной NVD (Normalized Variable Diagram). Для данной схемы гарантируется монотонность, однако ограничение по CBC вносит дополнительную диффузию. В [12] показано, что схема GAMMA слишком диссипативна для областей LES моделирования.
В настоящей работе предложена схема, также основанная на центральных разностях. Ее формулировка выполнена с использованием диаграмм NVD, однако, по сравнению с GAMMA, в ней применяется более "мягкое" ограничение, основанное на увеличении доли противопоточ-ности схемы лишь при увеличении отклонения нормализованной переменной от ее допустимого интервала. Также в схеме используется повышение уровня ограниченности в зоне пограничного слоя, позволяющее устранять численные осцилляции в зонах больших градиентов скоростей.
Для предложенной схемы показано, что на задаче конвективного переноса резкого фронта пассивного скаляра, она не приводит к возникновению численных осцилляций по сравнению с центрально-разностной, а также не "размывает" фронт как центрально-разностная схема с добавлением противопоточности. Сравнение численной диссипации предложенной схемы и калибровка для нее констант модели LES для свободных течений проведено на задаче о вырождении изотропной турбулентности. Калибровка для пристенных течений проведена на задаче о турбулентном течении в плоском канале, решение которой выполнено на последовательно сгущающихся расчетных сетках.
1. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ И ЧИСЛЕННЫЙ АЛГОРИТМ
Уравнения для расчета турбулентных течений методом LES могут быть получены из системы уравнений Навье—Стокса. Это осуществляется представлением переменных в виде суммы круп-
но- и мелкомасштабной составляющих и применением процедуры фильтрации (см. [2]). Для случая несжимаемой жидкости отфильтрованные уравнения имеют вид (см. [2]):
dui
д
дщ = 0 dx¡ _ dp
д
Р— + Р—u¡u¡ = -— + — (( + тsfs ). dt dx ¡ j dx¡ дхУ1 3>
(1)
В системе уравнений (1) используются следующие общепринятые обозначения: г, у = 1..3 — номер составляющей переменной по координатным осям, щ — скорость, р — плотность, р — давление, х,у — тензор вязких напряжений, который удовлетворяет реологическому закону Ньютона, устанавливающему связь между тензором вязких напряжений и тензором скоростей деформаций (см. [16])
мъ - 3 дг8 j
3 дхк
Sij =
дщ + dUj dxj dxi
где ц — молекулярная динамическая вязкость, 5 у — символ Кронекера. Тензор "подсеточных"
напряжений Т^ — дополнительный член, появляющийся при фильтрации уравнений. Его вычисление, как правило, базируется на использовании обобщенной гипотезы Буссинеска (см. [2]), аналогично тензору турбулентных напряжений при ЯАМВ подходе:
2 И sgs I Sij
1 dUk 5.
3 дхк 1
+ - к 5
^ 3*sgs ij'
Здесь ц , к— подсеточная вязкость и кинетическая энергия вихрей с размером менее размера фильтра А, в качестве которого обычно используется величина корня кубического из объема ячейки сетки (см. [2]). В качестве "подсеточной" модели турбулентности в работе рассматривается наиболее популярная модель подсеточной вязкости Смагоринского (см. [2]):
И sgs = Р (CsA)2 S, S = V 2SijSij,
где Cs — эмпирическая константа (константа Смагоринского).
(2)
Важная особенность метода LES — зависимость значения константы Cs (или иной, в случае использования другой подсеточной модели) напрямую от свойств используемой численной схемы, в частности от ее численной диссипации (см. [2], [9], [11]). Если диссипация численной схемы значительна, то значение константы Cs должно быть существенно уменьшено, поэтому, вообще говоря, для каждой численной схемы должна проводиться калибровка константы Cs. В случае пристеночных течений в выражение (2) вводится демпфирующий множитель, обеспечивающий нулевое значение подсеточной вязкости на твердых стенках (см. [2]):
Цsgs = Р(CsA)2 S
1 - exp
У_
v 25/
(3)
где y — безразмерное расстояние до твердой стенки.
Решение системы (1) осуществляется, как правило, с использован
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.