CC BY
156
А. Н. Нуриев, О. Н. Зайцева
РЕШЕНИЕ ЗАДАЧИ ОБ ОСЦИЛЛИРУЮЩЕМ ДВИЖЕНИИ ЦИЛИНДРА
В ВЯЗКОЙ ЖИДКОСТИ В ПАКЕТЕ OPENFOAM
Ключевые слова: осциллирующее движение, система уравнений Навье-Стокса, периодические режимы течения, метод
конечных объемов, пакет OpenFOAM.
В работе рассматривается моделирование задачи об осциллирующем движении кругового цилиндра в вязкой несжимаемой жидкости в пакете OpenFOAM. Большое внимание уделено вопросам построения численной модели в пакете: дискретизации области решения, выбору численной схемы, алгоритмов решения, оценке точности используемых методов. Представлены результаты расчетов в диапазоне низких чисел Стокса: построена карта режимов, выполнена визуализация характерных течений, проведено сравнение полученных данных с результатами лабораторных экспериментов.
Keywords: oscillatory motion, Navier-Stokes system of equations, periodical solutions, finite volume method, OpenFOAM package.
In this paper the modeling of the oscillatory motion of a circular cylinder in a viscous incompressible fluid in the package OpenFOAM is considered. Much attention is paid to the construction of the numerical model in the package: to discretization of the solution area, to the choice of the numerical scheme and algorithms of solution, to assessing the accuracy of the methods. The results of the calculations at the low Stokes numbers are presented: a map of regimes is constructed, a visualization of typical flows was done, a comparison of the data with the results of laboratory experiments was made.
Введение
В настоящее время прикладные вычислительные пакеты представляют собой мощный и хорошо развитый аппарат для решения широкого класса физических задач. Не секрет, что они активно используются ведущими мировыми инженерными компаниями при создании и испытании разнообразной высокотехнологичной продукции. Они позволяют экономить финансовые средства и время, заменяя множество дорогостоящих этапов разработки и тестирования численными экспериментами. Однако применение этих пакетов в области научных академических исследований до последнего времени ограничивалось двумя следующими факторами. Как известно, большинство вычислительных пакетов являются проприетарными программными продуктами и требуют крупных финансовых вложений на этапе их приобретения. Но более важно то, что они содержат закрытый программный код и как следствие представляют собой «черный ящик» для исследователя. Это ограничивает возможности по созданию, модификации и верификации новых численных моделей, мешает при оценке точности полученных результатов и т. д.
Появление на рынке свободного программного пакета ОреиРОАМ в последние несколько лет позволило изменить данную ситуацию. Широкий инструментарий для формализации задачи, высокая эффективность реализации, а также хорошая масштабируемость под архитектуру вычислительной системы позволяют легко сконструировать численную модель в пакете. Открытый исходный код в свою очередь дает возможность в деталях контролировать ход решения, начиная от построения сетки до выбора схем аппроксимации слагаемых управляющей системы и методов численного решения. Единственным ощутимым недостатком пакета
является его скудная документация, которая, несмотря на активное развитие пакета, пополняется в последнее время лишь статьями исследователей.
В настоящей работе рассматривается применение пакета ОреиРОАМ для решения задачи об осциллирующем движении кругового цилиндра в вязкой несжимаемой жидкости. В связи с указанным недостатком исследование сопровождается большим количеством информации о численной схеме, особенностях ее реализации в пакете и результирующей точности методов.
Задача об осциллирующем движении кругового цилиндра в вязкой несжимаемой жидкости относится к разряду перспективных областей исследования классической
гидромеханики. Задача хранит в себе большой нераскрытый теоретический и практический потенциал. Морское и гражданское строительство, авиационно-космическое проектирование,
робототехника - это лишь некоторые из областей, в которых задача имеет практическое приложение. С теоретической точки зрения огромный интерес представляет изучение сложных физических механизмов взаимодействия, структурных особенностей течения, анализ интегральных характеристик, исследование вопросов
устойчивости и бифуркаций решения. Еще один важный фактор, который привлекает современных исследователей к задаче - это обширная экспериментальная база результатов (например [1]-[3]), которая накопилась за несколько последних десятилетий. Она дает широкие возможности для верификации численных моделей и одновременно служит хорошей отправной точкой для разностороннего изучения задачи.
Различные подходы к численному
исследованию задачи ранее рассматривались в работах [4-8].
В безразмерной постановке задача
управляется двумя параметрами: числом Келигана-Карпентера КС и числом Стокса р. Первый характеризует отношение амплитуды колебаний к диаметру цилиндра, второй - квадрат отношения диаметра цилиндра к толщине нестационарного пограничного слоя. Параметры определяются соответственно:
и Т Б2
КС = _та^ , в =^_ .
Б УТ
Здесь итах- амплитуда скорости колебаний, Т -период колебаний, Б - диаметр цилиндра, V -кинематическая вязкость жидкости.
В данной работе исследования проводятся в диапазоне малых чисел Стокса в<55 при КС<10. Как показывают эксперименты, в этом диапазоне колебания цилиндра вызывают различные периодические режимы течения жидкости, структура которых зависит от обоих управляющих параметров.
Исследование двухпараметрической области, с большим количеством нестационарных решений требует огромного объема вычислений. Сложная структура режимов течения, неустойчивое поведение возле границ перехода, а также влияние трехмерных эффектов до сих пор не позволили провести полное моделирование задачи во всей области. В настоящей работе проводится большая серия расчетов (более 50 точек в параметрической плоскости) для плоскопараллельной модели задачи. Расчеты охватывают все зоны параметрической области, где экспериментально наблюдалось различное поведение течения.
Все вычисления проводятся на высокопроизводительном кластере в пакете ОрепРОАМ.
Конечными целями проводимого
исследования являются: моделирование
периодических режимов течения, локализация границ их устойчивости, а также верификация полученных результатов с экспериментальными данными.
В первом разделе работы рассматривается постановка задачи. Используемая схема дискретизации и алгоритмы решения представлены во втором, третьем, четвертом и пятом разделах работы. Шестой раздел посвящен вопросу моделирования возмущенного потока. В седьмом представлены основные результаты исследования.
Постановка задачи
Задача о движении осциллирующего цилиндра рассматривается в следующей постановке: круглый цилиндр с радиусом Я совершает высокочастотные гармонические колебания в вязкой несжимаемой жидкости по дифференциальному закону:
дх
колебания.
В рамках принятой модели процесс описывается нестационарной системой уравнений Навье-Стокса (1).
— + и -VU = -Vp + — AU dt Re , (1)
VU = 0
где U = (u, v) - безразмерная скорость, p -
безразмерное давление, Re = р КС - число Рейнольдса.
Для численного решения данной задачи удобно перейти в подвижную систему координат связанную с цилиндром. В этом случае для сохранения системы движения в форме (1) в новой неинерциальной системе координат, давление определяется как:
p = р + 2nfX sin(2rft)
Здесь первое слагаемое р - давление в подвижной системе координат, а второе вклад от инерциальных составляющих.
На границе цилиндра в новой системе координат задаются условия прилипания:
u =v = 0. (1.1)
На бесконечности изменение скорости определяются по следующему гармоническому закону:
u = cos (lrft\v = 0. (1.2)
В предположении о потенциальном течении жидкости на бесконечности из (1.2) можем получить также условия (1.3) для давления:
dp
дх
= 2nf sin(2f или p = 2nfxsin(lrft). (1.3)
—1 = - С05(2я/?), (о)
дґ
где - горизонтальная составляющая вектора перемещений, / = Кс - безразмерная частота
Дискретизация области решения
Для численного решения задачи рассматривается ограниченная область, которая представляет собой прямоугольный параллелепипед, в центре которого помещен круглый цилиндр (Рис
1.). Выбор трехмерной области для моделирования плоского течения проводится в соответствии с особенностью программного обеспечения. Размеры расчетной области указаны на Рис. 1. Граница области состоит из 7 частей: входная граница Г„ш и выходная граница Гои([е(, нижняя граница ГЪМот и верхняя граница Г(ор, передняя граница Ги задняя граница Гъаск , граница цилиндра . В
используемой декартовой системе координат ребра параллелепипеда параллельны основным осям, и плоскость течения параллельна плоскости хОу.
Для дискретизации расчетной области используются структурированные блочные сетки, построенные с помощью утилиты ЫоскМеБЬ, входящей в состав пакета ОрепРОАМ.
Область разбивается на непересекающиеся ячейки имеющие форму шестигранников. Разбиение в направлении оси 02 состоит из одной ячейки, так как в силу двухмерности задачи течение в этом направлении отсутствует. Также для учета
двухмерности на передней Г/Гоп( и задней ГЪаск
границах области, параллельных плоскости течения, используются специальные граничные условия, обсуждаемые в следующем разделе работы.
Рис. 1 - Расчетная область
В расчетах используются области двух размеров 50^30x1 и 50x60x1 с соответствующими сетками: т1, содержащей 105 ячеек и т2 — 2 -105 ячеек. Разрешающие особенности сеток
характеризуются следующим образом: объем самых мелких ячеек, расположенных на границе с цилиндром, составляет Ут = 0.004 для т1,
Ут = 0.001 для т2, количество ячеек на границе цилиндра равно Ыт = 80 для т1, Ыт = 160 для т2.
Ячейки наибольшего объема расположены в окрестности границ внешней области,
максимальный объем ячеек для т1 составляет
^тах = 0.02 , дЛЯ т2 - ^тах = °.°5 .
Помимо разрешающей способности важными факторами оценки качества построенных расчетных сеток являются [14]:
• неортогональность,
• скошенность,
• равномерность /х.
Эти характеристики оказывают ключевое влияние на точность результирующей аппроксимационной схемы.
Сетке т1 соответствуют следующие значения описанных характеристик: максимальная
неортогональность - 42.3828 градуса, средняя неортогональность - 4.11196 градуса. Максимальное значение коэффициента скошенности - 0.380392. Максимальное отклонение третьей характеристики /х = 0.4828 достигается в окрестности цилиндра вследствие сгущения в этой области расчетной сетки.
Сетке т2 соответствуют следующие значения характеристик: максимальная неортогональность -42.7689 градуса, средняя неортогональность - 3.1693 градуса. Максимальное отклонение /х = 0.476
также достигается в окрестности цилиндра.
Численная схема
Дискретизация системы уравнений движения (1) в пакете ОрепРОАМ проводится по методу
конечных объемов (РУМ) в декартовой системе координат. Для этого дискретные значения составляющих скорости и дискретные давления локализуются в центрах ячеек построенной расчетной сетки.
Для произвольной ячейки сетки с объемом У система уравнений (1) записывается в следующей интегральной форме:
(2)
V•UdV = 0
Первое слагаемое системы аппроксимируется в центре ячейки как произведение среднего значения подынтегральной функции на объем ячейки V. Для вычисления остальных объемных интегралов по контрольному объему V системы уравнений (2) используется общая процедура Гаусса [9], согласно которой осуществлялся переход от объемного интеграла к поверхностному. Далее поверхностные интегралы представляются в виде суммы интегралов по граням ячейки и приближенно вычисляются по формуле средних прямоугольников. После этого, полудискретная система уравнений для произвольной ячейки представляется в виде:
^1 +1 -{ии)/ =-! Я,р,
д(
/
+ у£^ •{Vи)
/
(3)
'/’
у Я/-{и) / = 0.
//
/
Здесь индекс / указывает на то, что переменная или градиент определены на грани ячейки, Р - в центре ячейки, а Я/ определяется как вектор ортогональный к грани ячейки и по модулю равный площади этой грани.
Для линеаризации системы (3), конвективные слагаемые представляются в следующей форме:
у Я/ •{ии)=у р{и), // где Р - массовый поток через грань с индексом / считается известным. Обновление Р связано с итерационной процедурой метода решения задачи.
Значения функции и нормальных градиентов на поверхности ячеек в системе (3) для внутренних ячеек расчетной области интерполируются из значений функции в центрах соседних ячеек. Рассмотрим применяемые в данной работе схемы интерполяции переменных.
Для аппроксимации градиента давления в расчетах применяется линейная интерполяция. Значение р/ на грани / между двумя ячейками с
центрами Р и N (рис. 1.) находится по формуле:
Р/ = Ррух + PN {1 - Ух К где / х - интерполяционный фактор, определенный выше. Порядок точности используемой аппроксимации обусловлен локальными
р
характеристиками сетки. Он понижается до первого в случае локальной скошенности сетки, т. к. интерполируемое значение определяется не в центре грани. В остальных случаях аппроксимация имеет второй порядок точности.
Для интерполяции переменных в конвективных слагаемых используется нелинейная NVD (normalised variable diagram) схема «Gamma», предложенная в работах [12, 13]. Применение этой схемы позволяет обеспечить устойчивость численной задачи, внося минимальную численную диффузию. В качестве схемы высокого порядка точности в «Gamma» используются линейная интерполяция, в качестве безусловно устойчивой схемы низкого порядка - противопоточная схема «upwind». Вычисление дискретных составляющих скорости (uf, vf) = Uf на грани с индексом f для F >
0 производится по формулам:
uf = uP
uP < О, uP > 1
Pm < uP < 1
V = uPfx + uN (l - fx ) uf =(1 -ru (l - fx ))uP +ru I1 - fx К , 0 < ~p <pm
vf = vp
vP < 0, vP > 1
vf = vPfx + vN (l - fx ), Pm < ~P < 1
Lvf =(1 -rv(l - fx ))vP +rv(1 - fx )vN ,0 < ~P <pm
Здесь иР =1 -((Уи)/ •d)/(2{'^7и)Р •dК,
~Р = 1 -(СУ^у)/ • d)/(2(Vv)р • d) - нормализованные
переменные, d - вектор направленный из точки Р в точку N Г и = ~Р / рт , Гу = ~Р / рт - факторы смешивания, 1/10 <Рт < 1/2 - предопределенная константа метода. Выбор больших значений Рт из указанного диапазона обеспечивает наилучшую устойчивость схемы, меньших - увеличивает точность интерполяции. В данной работе использовалось значение рт = 0.25 . В случае противоположного направления потока Р Р<0) формулы изменяются соответствующим образом.
В диффузионных слагаемых при дискретизации оператора Лапласа необходимо вычислять нормальные градиенты скорости на поверхности ячейки. На ортогональных участках сетки, где вектор Я параллелен вектору d, они вычисляются из значений скорости в центрах соседних ячеек по симметричной схеме второго порядка:
S .(Vu )f = |4
|d|
(4)
На неортогональных участках сетки скалярное произведение £ .(Ум) / представляется в виде суммы двух слагаемых:
£. (у и)/ = I • (Ум )/ + к • (Ум )/ .
Первое из слагаемых отвечает за ортогональный вклад, второе - за неортогональную поправку, при этом для векторов к и I выполняется соотношение
£ = к +1. (5)
Ортогональный вклад вычисляется по формуле (4), где вместо вектора £ используется вектор I, параллельный вектору С, длинна которого определяется по формуле:
l=
d
d • S
IS2.
Неортогональная поправка вычисляется следующим образом: вектор к находится из соотношения (5), а
V) /
значение
ячейки
градиента / на грани
интерполируется из значений градиентов в центрах соседних ячеек:
(Ум )/ = (Ум )р/х +(м )„ (1 - /х ).
Для интерполяции слагаемого (уи)/ применяется
аналогичный подход.
Аппроксимация диффузионных слагаемых имеет второй порядок точности для равномерных участков сеток, на неравномерных участках порядок понижается до первого [14].
Для дискретизации системы (3) по времени используется неявная схема Эйлера:
V+ZFun=-У Sfpf ZSf .(VU )n/
(6)
+v
f
u О=0
(7)
У£/р}=0 У£/•
/ ’ /
Здесь верхние индексы «о» и «}» указывают на использование переменной со старого или нового временного слоя соответственно, Т - шаг по времени. Хотя сама схема является безусловно устойчивой, но для минимизации эффектов связанных с аппроксимацией первого порядка точности, шаг по времени во всех расчетах выбирается из условия С™* < 0.1 - максимальное число Куранта не превышает значения 0.1. Число Куранта в пакете ОрепРОЛМ определяется по формуле:
= Рр \Т
Co =
S
где |иР| - модуль скорости в ячейке, 8 - размер ячейки в направлении вектора скорости.
Граничные условия
Для замыкания результирующей дискретной системы уравнений необходимо определить граничные условия для расчетной области (рис. 1.) и провести дискретизацию в граничных ячейках.
На входной Г1пШ и выходной ГоШ1е1 границах области задаются неотражающие граничные условия вида:
р = 2ж/х б1п(2^/?), и0 > 0
— = 2fx sin(2ft), uQ < 0
дх
v=0
Они скомбинированы из условий (1.2), (1.3),
определенных на бесконечно удаленной границе. Величина горизонтальной составляющей скорости вычисляется при решении задачи. Условие для
т
N P
давления зависит от переменной u0 , определяющей направление потока относительно внешней нормали к границе. Она равна u0 = cos(2^/t) для входной границы, uo = - cos(2^/t) - для выходной границы.
На верхней rtop и нижней Fbottom границах ставятся условия проскальзывания:
dp „ du п
-*- = 0, — = 0, v = 0.
dy dy
Эти условия также являются следствием условий (1.2), (1.3) на бесконечности.
На границе цилиндра Гcylinder ставятся условия прилипания для скорости: u=v=0, и условие для давления:
^ = 0.
dn
На передней Гfront и задней rback границах
области задаются специальные «пустые» граничные условия, предусмотренные в пакете для случаев, когда вычисления в обозначенном направлении не проводятся.
Все представленные граничные условия делятся на два вида: условия Дерихле и условия Неймана. Дискретизация уравнений движения (3) во всех граничных ячейках сводится к рассмотрению этих двух различных случаев.
При аппроксимации конвективных слагаемых на границе с условиями Дирихле определены значения скоростей ub , vb , что без дополнительной интерполяции позволяет вычислить
соответствующие граничные компоненты:
Fub , Fvb .
Аналогично для аппроксимации градиента давления в случае граничных условий Дирихле используется значение pb на грани ячейки.
Для аппроксимации диффузионных слагаемых на границах с условиями Дерихле необходимо вычислять нормальные градиенты скорости. Поскольку используемые сетки всюду ортогональны границам расчетной области, т. е. вектор d, направленный из центра любой граничной ячейки в центр ее грани, принадлежащей границе области, параллелен вектору нормали к этой грани, соответствующие нормальные градиенты могут быть вычислены по формулам:
ub — up
S .(Vv) = 14
'b VP
и щ ■
Для случая граничных условий Неймана значения скоростей иъ, уъ в центре грани ячейки, необходимые для аппроксимации конвективных слагаемых, могут быть вычислены по следующим формулам:
ub = uP +
d .(Vu)
b , vb = vP +
d. (Vv )b
где и -{Уи )ъ, и -(у) - известные нормальные градиенты.
По аналогичной формуле, в случае граничных условий Неймана, вычисляется значение
pb на грани ячейки для аппроксимации градиента давления:
pb = pp + d • (vp)b.
При аппроксимации диффузионных
слагаемых на границах с условиями Неймана значения нормальных градиентов скорости известны:
£ -(Vu )b, £ • (Vv)b .
Для реализации чередующихся граничных условий на входной rinlet и выходной routlet границах используется дополнение groovyBC [17]. В качестве начальных условий задачи во всей расчетной области используются значения скоростей и давления соответствующие невозмущенному
потоку.
Численное решение
Дальнейшее решение дискретной задачи в пакете основано на подходе («segregated approach») раздельного решения уравнений для скорости и давления. Для соблюдения согласованности аппроксимационных схем, уравнение для давления выражается из дискретизованных уравнений движения и неразрывности (6) - (7).
Представим уравнение движения в
следующем полудискретном виде:
apUP =-У CnUN + Up -Vp .
(8)
Недискретизованным здесь остается градиент давления, который определен как произведение среднего значения подынтегральной функции на объем ячейки и локализован в центре ячейки. Дискретизация других слагаемых проведена согласно представленным выше схемам (табл. 1), коэффициенты при соответствующих компонентах
скоростей Рр и РN содержатся в диагональных матрицах ар и ан размерности 2x2. В приведенной форме (8), все слагаемые уравнения разделены на объем ячейки V.
Выразим из уравнения (8) вектор скорости в центре расчетной ячейки ир :
ир = арИ(и)-арУр, (8.1)
где за И(Р) обозначен оператор вида:
И (и )=-у аыипы + .
N
Полученный вектор скорости можно
интерполировать на грань ячейки следующим образом:
Р/ = (а^ И (и))/ - (ар1)/ (Ур)/ , (9)
где (ар*)/ - обозначает интерполяцию
соответствующих коэффициентов матрицы.
Подставляя (9) в уравнение неразрывности, получим дискретное уравнение для давления:
У £/ .{ар)/(Ур)/ =у £/ .{арИ(и)
//
Аппроксимация оператора Лапласа здесь
N
проводится так же, как и в диффузионном члене.
Результирующую систему уравнений с выделенным уравнением для давления запишем в виде:
ари р = И {и) - у £/р/
коррекции потоков:
f
(10)
У sj
f
( 1 ^ — Vp
(
H (U )
Л
/f
(11)
(11)
= Е Sr
7 f f
Для решения задачи в форме (10) используется программа icoFoam, реализующая алгоритм PISO (Pressure Implicit Splitt Operator) [9] -[10]. Итерационная процедура алгоритма основана на последовательном решении уравнений для скоростей и давления. Согласно используемой реализации [13] вычисление неизвестных полей на новом временном слое проводится по следующей схеме:
0. По полю скоростей U, известному с предыдущего временного слоя, вычисляются потоки F.
1. Проводится неявное вычисление предиктора нового поля скоростей. Для этого решается система, образованная из уравнений (12), где поле давления берется со старого временного слоя:
apUP + £ aNUnN = Up -£ Sfp}
N Т f
(В общем случае найденный предиктор не удовлетворяет уравнению неразрывности).
2. Проводится k-циклов коррекции:
2.1. По найденному приближению поля скоростей вычисляются значения оператора H(U):
UP
H (U )=-У CnUN +
N
*
и рассчитываются значения и Р : и*Р = аР1И (и ), представляющие собой новые значения скоростей без учета давления (см. (8.1)).
2.2. Найденные скорости интерполируются на границы и*Р ^ и**, после чего вычисляются соответствующие значения потоков:
F = Sf —f +
(fo-sj -P p (■),.
Второе слагаемое в правой части выражения
служит для замены слагаемого
Sf. uf (-
(cp1 ff
оператора
(cp1H (U ))
f, полученного в результате
F<o / \
интерполяции, слагаемым ---------(1)
т
- известным с
предыдущего временного слоя, для которого в точности выполняется уравнение неразрывности.
2.3. Рассчитывается новое поле давления. Для этого решается система уравнений вида:
У sj -(ap1 f(VpPj =У
Fff.
ff 2.4. Найденные давления используются для
F = F * + Sf
.(cp1 ff (vp)j
>/ Гр }/\УН}/
и поля скоростей:
и„ = и*Р + а-} (Ур)
2.5. Выполняется коррекция граничных условий.
В представленном алгоритме следует отметить несколько важных моментов:
• При малых числах Куранта связи между скоростью и давлением в уравнениях движения считаются более сильными, чем нелинейные связи на фиксированном временном слое. В связи с этим матрицы коэффициентов аР, ам^, зависящие от значений потоков Р, обновляются на каждой итерации по времени только на этапе предиктора.
• В алгоритме решаются две линейные задачи: первая для вычисления предиктора поля скоростей (шаг 1), вторая для расчета нового поля давления в цикле коррекций (шаг 2.3). Все остальные операции выполняются по явным формулам.
• При решении линейных задач для обеспечения диагонального преобладания матриц соответствующих систем используется отложенная коррекция [15]. Согласно этому методу конвективные слагаемые представляются в виде суммы:
(
У F—)f =У F (u )}p
ff
+
У
v f
F (- )f-У
f
F (- )
Первое слагаемое в правой части, полученное в результате дискретизации с использованием противопоточной схемы «upwind», не нарушает диагонального преобладания матрицы и поэтому используется для определения коэффициентов системы. Оставшийся комплекс, представляющий поправку более высокого прядка, определяется явным образом и добавляется в столбец свободных членов.
Аналогичный подход используется для дискретизации диффузионных членов на неортогональных участках сетки. Как было отмечено ранее, дискретизация этих слагаемых представляется в виде:
Е Sf • (Vu )f = £ l (Vu )f + £ k • (Vu )f f f f
Ортогональный вклад используется для определения коэффициентов системы,
неортогональная поправка задается явным образом и добавляется в столбец свободных членов.
Для решения линейных систем с необходимой точностью, с учетом вышеописанного
представления слагаемых, требуется несколько итераций: решение, полученное на каждой
итерации, используется для обновления добавочных слагаемых в столбце свободных членов. В случае системы уравнений для давления такой цикл коррекций реализуется на шаге 2.3. Точное решение задачи здесь обеспечивает консервативность результирующих потоков.
Основные параметры алгоритма PISO - число коррекций к (шаг 2) и число неортогональных
т
коррекций (шаг 2.3) - выбираются равными
соответственно (3, 3) для используемых сеток. Для решения системы уравнений для давления (шаг 2.3) применяется метод сопряженных градиентов PCG c геометро-алгебраическим многосеточным
предобуславливателем GAMG (Geometric agglomerated algebraic multigrid solver). В реализации GAMG для сглаживания используется метод Гауса-Зейделя с числом 1, 2 пре- и пострелаксаций соответственно, для агломерации ячеек сетки используется faceAreaPair алгоритм [16]. Система уравнений для скоростей решается методом би-сопряженных градиентов PBiCG с предобуславливателем основанным на неполной LU факторизации (Diagonal incomplite-LU). Сходимость по всем методам выполняется до значений невязки меньших 1e-8.
Все вычисления выполняются
распределенным образом по технологии MPI, применяя метод декомпозиции области решения (domain decomposition). Для этого расчетная область делится на 2-4 подобласти по вертикали в зависимости от рассматриваемого случая. Подзадачи в каждой подобласти рассчитываются на различных ядрах процессора.
Моделирование возмущенного потока
До некоторого критического значения KCcr при фиксированном значении параметра в течение в области остается устойчивым к внешним возмущениям и сохраняет симметрию и периодичность. При значениях KC больших критического устойчивость теряться, и течение утрачивает исходную симметрию под действием внешних возмущений. В случае реального физического явления или при проведении лабораторных экспериментов такие возмущения всегда присутствуют. В случае прямого численного моделирования природа естественных внешних возмущений может быть связанна только с вычислительными ошибками или ошибками
дискретизации. При проведении численных экспериментов в области умеренных КС > КСcr влияние вышеописанных ошибок на симметрию проявляется не раньше, чем после 15-20 периодов колебания.
Для моделирования переходных режимов в окрестности KCcr естественных возмущений оказывается мало, поэтому в область течения вносятся дополнительные возмущения по методу Мартинеза [18]. Этот метод ранее применялся
Джастинсоном для решения задачи о гармонических колебаниях цилиндра в работе [6]. Согласно этому методу возмущения порождаются вращением
цилиндра относительно оси Oz: некоторый
промежуток времени цилиндр поворачивается по часовой стрелке, а затем, через небольшой временной интервал, он поворачивается в противоположном направлении до исходной позиции.
Результаты
Все проведенные расчеты отмечены на карте
режимов (рис. 2) в параметрической плоскости в-КС. В зависимости от наблюдаемой структуры течения они обозначены различными маркерами.
В большей части диапазона преобладают периодические режимы. По характерным структурным различиям можно выделить шесть типов периодического течения. Ромбовидными и круглыми маркерами обозначены симметричные относительно плоскости колебания режимы. Они реализуются в области малых чисел КС. Характерная структура течения для симметричных режимов, представленная на рисунке 3 (верхний). Для визуализации использовались невесомые частицы, которые вымывались течением из малой окрестности цилиндра.
КС
10
8 6 4 2 0
Рис. 2
При переходе через КСсг (Р) симметричные режимы теряют устойчивость. Потеря устойчивости происходит в пользу несимметричных периодических режимов, отмеченных на карте квадратными и крестообразными маркерами. Течение отклоняется от оси колебания цилиндра. В случае режима, обозначенного крестообразными маркерами, существенно изменяется и период течения. Зоны устойчивости этих переходных режимов располагаются вдоль границы потери симметрии. Характерные картины течения представлены на рисунке 3 (нижний) и рисунке 4 (верхний).
Рис. 3 - Структура течения для р=35, KC=4.5 и в=35, KC=5
С дальнейшим ростом КС течение переходит в псевдопериодический режим, обозначенный треугольным маркером (▼). В этой зоне слабая
А А \ +
\ \
\. F \ + +
Xd "’V, ▼ 7.......
° ▼ Е
__□ ▼
X
О——i “ —
О ----5
о
о С В о о о -
t А* 1 о о
15 20 25 30 35 40 45 50 (3
■ Карта режимов
периодичность наблюдается в течение нескольких циклов, затем происходит смена направления сброса вихрей. С увеличением КС периодичность полностью утрачивается. Точки с непериодическим режимом течения обозначены на карте маркерами +.
В области больших КС был локализован еще один периодический режим течения (отмечен маркером ▲). Режим характеризуется формированием диагональных вихревых дорожек. Они образуются из вихревых пар, которые каждые полпериода срываются с разных сторон цилиндра (рис. 4 нижний).
Рис. 4 - Структура течения для р=35, KC=5.5 и в=35, KC=8
Полученные данные хорошо согласуются с экспериментальными наблюдениями. Аналогичные результаты в частности были получены в экспериментальной работе Татсуно и Бирмана [1]. На рис. 1 пунктирными линиями отмечены зоны периодических режимов полученные этими авторами. Можно установить следующее соответствие между обозначениями: А* - ◊, А - о, С -X, Б - □, Е - ▼, Б - ▲.
Как видно большинство моделируемых точек лежат в границах экспериментальных режимов.
Хорошая согласованность результатов наблюдается и в области перехода между симметричными и несимметричными режимами.
Небольшое смещение зоны режимов наблюдается в окрестности границы режима B. Для плоско-параллельной задачи режим B с преобладающими, согласно Татсуно и Бирману, трехмерными свойствами становится неразличимым. Как видно на рисунке его место, согласно представленной двухмерной модели, занимает режим A. Что касается области режима G, схема которого представлена в работе [1], то его место занимают непериодические течения.
Литература
1.M. Tatsuno et al., J. FluidMech, 211, 157-182 (1990).
2. H. K. Williamson, J. Fluid Mech., 155, 141-174 (1985).
3. T. Sarpkaya, J. Fluid Mech., 165, 61-71 (1986).
4. H. Dutsch et al., J. Fluid Mech., 360 249 - 271 (1998).
5. B. Uzunoglu et al., International Journal for Numerical Methods in Engineering, 50, 2317-2338 (2001).
6. P. Justesen, J. Fluid Mech., 222, 157-196 (1991).
7. А.Н. Нуриев, Вестн. Каз. технологического ун-та., 16, 334-336 (2011).
8. А.Н. Нуриев и др., Вестн. Каз. технологического ун-та., 4, 104-109 (2013).
9. J.H. Ferziger et al. Computational methods for fluid dynamics, Springer, Berlin, 2002, 424 p.
10. H.K. Versteeg et al., An introduction to computational fluid dynamics. The finite volume method, Longman, New York, 1995, 258 p.
11. OpenFOAMs basic solvers for linear systems of
equations. URL: http://www.tfd.chalmers.se/~hani
/kurser/OS_CFD_2008/TimBehrens/tibeh-report- fin.pdf
12. H. Jasak et al., Int. J. Numer. Meth. Fluids, 31, 431-449 (1999).
13. H. Jasak, PhD. Thesis, Imperial College, University of London, 1996, 394 p.
14. Franjo Juretic, Error, PhD. Thesis, Imperial College, University of London, 2004, 233.
15. P. К Khosla et al., Computers & Fluids, 2, 2, 207-209 (1974).
16. Open FOAM (The Open Source CFD Toolbox): User Guide Version 2.1.1, 2012 (http://www.openfoam.org/docs/).
17. URL: openfoamwiki.net/index.php/Contrib_ groovyBC
18. G Martinez, These docteur-Ingenieur, I.N.P. Toulouse., 1979.
© А. Н. Нуриев - асс. каф. информатики и прикладной математики КНИТУ, асп. каф. аэрогидромеханики КФУ, 2anuriev@kpfu.ru; О. Н. Зайцева - асс. каф. информатики и прикладной математики КНИТУ, olga_fdpi@mail.ru.
