CC BY
44
Вычислительные технологии
Том 16, № 5, 2011
Ускорение сходимости решения нестационарных задач динамики несжимаемой жидкости*
И.Ф. Ешкунова, С. Г. Чёрный, Д. В. Чирков Институт вычислительных технологий СО РАН, Новосибирск, Россия е-таП: cher@ict.nsc.ru, dchirkov@ngs.ru
Рассматриваются различные подходы к ускорению сходимости решения нестационарных задач динамики несжимаемой жидкости методом, предложенным авторами ранее. В его основу положены концепции искусственной сжимаемости, конечных объемов и приближенной факторизации. Сравнение представленных подходов проводится на задаче моделирования течений в проточном тракте реальной гидротурбины.
Ключевые слова: нестационарные течения, несжимаемая жидкость, численное моделирование, ускорение сходимости, метод искусственной сжимаемости, метод приближенной Ы1-факторизации, гидротурбина, прецессия вихревого жгута.
Введение
В настоящее время для решения трехмерных уравнений несжимаемой жидкости Эйлера, Навье — Стокса и Рейнольдса, записанных в простейших переменных, используются два класса методов: расщепления, или проекционные, и искусственной сжимаемости. В первых используется идея о возможности определения давления на последующем временном слое из условия обращения в нуль дивергенции вектора скорости [1-4]. При этом возникает необходимость решения на каждой итерации текущего момента времени уравнения Пуассона для давления. Время, затрачиваемое на его решение, составляет значительную часть от общего времени выполнения итерации.
Второй класс методов базируется на концепции искусственной сжимаемости [5-8], заключающейся во введении в уравнение неразрывности производной от давления по псевдовремени, а в уравнения количества движения — производных также по псевдовремени от соответствующих компонент скорости. На каждом шаге по времени организуется процесс установления по псевдовремени.
Метод искусственной сжимаемости (МПС) взят за основу предложенного в [9] метода решения трехмерных уравнений несжимаемой жидкости. В нем на каждой итерации по псевдовремени необходимо решать систему линейных алгебраических уравнений (СЛАУ), возникающую в процессе дискретизации основных уравнений. Таким образом, экономичность построенного в [9] метода решения уравнений несжимаемой жидкости напрямую зависит от (а) затрат на решение СЛАУ, (б) количества итераций по псевдовремени.
В работе [9] для уменьшения вычислительных затрат на этапе а) исходная СЛАУ приближенно подменяется факторизованной, в которой матрица представляется в ви-
* Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00475-а).
де произведения двух экономично обращаемых треугольных матриц: нижней и верхней, Для реализации краевых условий за каждой границей расчетной области вводятся два ряда фиктивных ячеек. Такой подход позволяет использовать общие процедуры вычисления разностных потоков для всех внутренних ячеек расчетной области. Величины искомых параметров в фиктивных ячейках в СЛАУ не входят и заполняются в соответствии с типом краевых условий после вычисления их значений внутри области.
При обращении треугольных матриц невязки в фиктивных ячейках полагаются нулевыми, В силу этого будем говорить, что в [9] используется явная процедура реализации граничных условий, В результате достигнуто весьма экономичное обращение факторизованной матрицы в СЛАУ, но не изучены следующие вопросы:
— как приближенная 1Д -факторизация влияет на время установления по псевдовремени;
— уменьшится ли количество итераций, если, разрешая СЛАУ, более точно аппроксимировать исходную матрицу в ней;
— как повлияет на сходимость итераций неявная реализация краевых условий,
В настоящей работе изучение сформулированных вопросов проводится на задачах моделирования пространственных течений в проточном тракте (ПТ) гидротурбины (ГТ). В качестве тестовых задач рассмотрены стационарное и нестационарное течения в конусе отсасывающей трубы (ОТ), Приведены сравнения требуемых итераций и времен, затраченных на расчет.
Еще один способ ускорения расчетов состоит в применении блочного распараллеливания, Расчетная область разбивается на несколько блоков, в каждом из которых расчеты проводятся одновременно на отдельных процессорах, после чего между смежными блоками происходит обмен граничными значениями. Приводятся данные ускорения расчета течения в проточном тракте ГТ, Результаты расчетов сравниваются с экспериментальными данными.
1. Основные уравнения и постановки задач
1.1. Уравнения Эйлера и Рейнольдса движения несжимаемой жидкости
В декартовой системе координат (х^х2,х3) уравнения Эйлера и Рейнольдса можно записать в следующем обобщенном виде:
дм,
дх,
0,
д / ч д
дх,
д
-р
дх
дм, дм,
Р
р
р
3
2-к зк'
дх, дх,
+ ¡г
(1) (2)
где и = (м1,м2, мз) — декартовы компоненты вектора скорости; р - гидростатическое давление; р — плотность; к — кинетическая энергия турбулентности; иед — коэффициент эффективной вязкости; ¡г — массовые силы. Например, если вращение системы координат с постоянной скоростью ш происходит вокруг оси Ох3, а сила тяжести действует в положительном направлении той же оси, то ¡1 = х1ш2 + 2м2ш, ¡2 = х2ш2 — 2м1ш, ¡3 = д, где д — ускорение свободного падения.
Уравнения Эйлера получаются из (1), (2) при к = 0 и г/ея = 0,
1.1.1. Векторная форма записи основных уравнений
Для удобства дальнейшей работы с основными уравнениями запишем их в виде единого векторного уравнения в дифференциальной форме
+ Е Х1
где
СтЖ2 + НХ3 ,
<2 = (Р,«1,И2,из)Т, к4 = ^аё(0,1,1,1), Е = (0,/1,/2,/э)Т,
/
Е
и1
\
и 1 + р - Т11 «1«2 - Т21 \ и1из - Т31 )
(
С
и2
\
Т12 Т22
\ П2Щ - Т32 )
и2и1 П22 + р
Н
/ из
изи1 - Т13 ПзП2 - Т23
\ из2 + р - Тзз )
ди,- ди,
Ч]
дх ] дх,
(3)
1.1.2. Интегральная форма
Для построения метода конечных объемов в работе используется эквивалентная дифференциальному уравнению (3) на гладких решениях интегральная форма этого уравнения
Ы4
д_
т
<3 ¿V = - Л к^э
I е ¿у (4)
V дУ V
где дУ — замкнутая поверхность произвольного фиксированного объема V, ^ = п •
^Б = ^Б2, ^Бз) — элемент поверхноети Б, умноженный на единичную внешнюю нормаль п к ней. Вектор потоков представляется как сумма невязкого и вязкого потоков
К • ^ = (К,п + К
(5)
где
,
/ и \
^и + ре 1 и^и + ре 2 \ изи + рез )
Кг
000
Т11 Т12 Т13 Т21 Т22 Т23 \ Т31 Т32 Т33 )
е1 = (1,0,0), е2 = (0,1,0), ез = (0,0,1) - базис декартовой системы координат.
1.2. Модель турбулентности для уравнений Рейнольдса
Для замыкания уравнений Рейнольдса используется к - е-модель турбулентности Кима—Чена [10]. Турбулентная кинетическая энергия находится из уравнения переноса к
дк д Л дк \ _ 2^к
д£ дх,-
ки, - ик
дх
С - е -
],
У2
е
де д д£ дх,
де
I дх,
е е2
= Се\ — С — Се 2 — ; кк
(6)
где и£ = V + (Те
1.15 С£1 = 1.15 С
£2
1.9.
1.3. Постановки задачи численного моделирования течения в гидротурбине
При численном моделировании течения в ГТ ввиду сложности ее геометрии и наличия неподвижных проточных частей (спираяьпая камера, каскад статорпых колопп с лопатками направляющего аппарата (НА), ОТ) и вращающихся вместе с рабочим колесом (РК) межлопастпых каналов проточный тракт необходимо разбивать па отдельные блоки (рис. 1).
Моделирование течения в ГТ можно осуществлять в различных приближениях: в полной и циклической постановках (см. рис. 1), во всем проточном тракте ГТ или в некоторых его сегментах [9], а также в стационарной и нестационарной постановках. Наиболее экономичной с точки зрения вычислительных затрат является циклическая постановка, в которой принимается допущение, что течения во всех межлопаточпых каналах НА и РК циклически идентичны. Расчет проводится только в одном из каналов НА и РК, а па боковых границах каналов ставятся условия периодичности потока. Для передачи параметров потока из вращающихся сегментов в неподвижные и наоборот их значения усредняются в окружном направлении.
1.4. Краевые условия
Для реализации краевых условий за каждой границей всех блоков расчетной области вводятся два ряда фиктивных ячеек, значения искомых переменных в которых па каждой итерации метода установления заполняются в соответствии с типом границы. Такой подход позволяет использовать общие процедуры вычисления разностных потоков для всех внутренних ячеек расчетной области.
1.4.1. Входная и выходная границы
Во входном сечепии расчетной области задается распределение скорости и, в выходном — давления р, а также, в случае вязких течений, касательная компонента скорости.
Рис.. 1. Сегментированы!! проточный тракт гидротурбины: а полная постановка; б цикли ческая постановка
1.4.2. Твердая стенка
На твердых стенках в случае невязких течений ставится условие непротекания и • п|ш = 0, в оду чае вязких течений — условие п рилипания и|ш = и где — скорость движения стенки.
1.4.3. Граница между соседними блоками
В случае течения в области, состоящей из нескольких блоков, расчет проводится следующим образом, В каждом блоке последовательно осуществляется по одной итерации метода установления, после чего проводится обмен всеми параметрами течения на той части их общей границы, через которую протекает жидкость. Затем расчетные итерации повторяются. Сетки строятся таким образом, чтобы можно было легко осуществлять интерполяцию данных из рабочих ячеек в фиктивные между смежными блоками.
1.4.4. Граница периодически повторяющегося течения
В циклической постановке, в которой расчет проводится только в одном из каналов НА или РК, на боковых границах каналов ставятся условия периодичности потока. Значения параметров потока в других межлопаточных каналах НА и РК могут быть получены с помощью нескольких поворотов вычисленных векторов скорости и давления из ячеек расчетной области вокруг оси РК на угол периода. Угол периода определяется количеством лопаток в НА или лопастей в РК,
2. Численный метод
2.1. Метод искусственной сжимаемости
Уравнения (3) решаются численно с использованием подхода искусственной сжимаемости и неявного метода конечных объемов [9]. Следуя идее метода искусственной сжимаемости [5, 6], в уравнение неразрывности добавляется производная от давления по псевдовремени
дР | рдщ
дТ дх,
0,
(8)
где в _ коэффициент искусственной сжимаемости. Производные от компонент скорости также добавляются в уравнения сохранения момента импульса
д др д <9х,- 3 дхдxi
' дщ ди/ дх, дх,
+ /,.
(9)
Модифицированная система уравнений (8), (9) решается численно. При этом на каждом шаге по физическому времени £ проводится установление по псевдовремени т.
Данный подход может быть использован для решения как стационарных задач методом установления, так и нестационарных.
2.2. Неявная конечно-объемная аппроксимация
Для построения консервативной разностной схемы модифицированные уравнения (8), (9) представляются аналогично (4) в форме интегральных законов сохранения для произвольного фиксированного объема V:
+ К'-^Л I ЦёУ + у Н • Ж = I ТёУ, (10)
^ ' V дУ V
где ^ — поверхность объема V, ¿Б = п • ¿Б — элемент поверхности, умноженный на единичный нормальный вектор, Q = (р, щ, и2, из)т, Е = (0, /1, /2, /з)т,
( ви1 ви2 ви3 \
Н =
и? + р — Т11 ЩЩ — Т12 и1из — Т13
2 у-
и1и2 — Т12 и? + Р — Т22 и?из — Т23
2
\ и1из — Т13 и2из — Т23 и2 + Р — Т33 /
Ит = ^(1,1,1,1), И* = ё1а§(0,1,1,1).
Введем следующие обозначения для векторов О и Е, осредненных по ячейке г^'к с объемом Vijk\
щк = ! щк = ! ¡гм,
У1зк У1зк
и отнесем их к центру ячейки. Здесь верхний индекс п — номер слоя по времени. Дискретизацией уравнения (10) получим
- (сп+1)^+1 — (Сп+1)^ з(Сп+1)5+1 — 4Оп + 1
К -:--Ь К-:-
А г 2А£
V = (ИНБ^1)^1, (11)
где Дт и Д£ — шаг то псевдовремени и шаг по времени соответственно, 5 — номер итерации по псевдовремени. Правая часть (11)
ИНБ = — (К • Б)т/2 — (К • S)i—1/2 + (К • Б),.+1/2 — (К • Б)^-1/2 + + (К • Б)к+1/2 — (К • Б)к—1/^ + Е^ где (К • S)i+l/2, (К • (К • Б)
к+1/2 разностные тотоки через грани г + г^' + 1/2к и г^'к + 1/2 ячейки с номером г^'к и объемом Vijk■
Разностный поток представим в виде суммы конвективного (невязкого) и вязкого потоков. Невязкие разностные потоки вычисляются таким образом, чтобы результирующая схема имела направленные против потока разности второго или третьего порядка аппроксимации. Для обеспечения монотонности решения в случае невязких течений или при больших числах Рейнольдса в случае вязких течений применяются ограничители невязких потоков, обеспечивающие выполнение для разностной схемы свойства ТУБ, Аппроксимация невязкого потока Кт^ • Бт+1/2 через грань Бт+1/2 ш-й ячейки сетки
осуществляется в соответствии со схемой Чакраварти — Ошера [7]: к*™.8 1
т+1/2
К™ (От) + К™ (Ст+1) Бт+1/2 — |Ат+1/2|Дт+1/2С — т+1/2, (12)
д (К™ (О С)
где А(С*) = 1 ^ ' А = А+ + А", |А| = А+ - А", Ат+1/20, = С*т+1 - <Эт, а Wm+1/2 обеспечивает третий порядок аппроксимации схемы.
2.3. Линеаризация
Уравнение (11) линеаризуется с использованием метода Ньютона
(QS+1 - Qs) =
(—Rr + Л-В*) V - (¿RHS
VAr 2At J \dQ
= -r,3Q'~42X+q""+rhs'- («)
здесь верхний индекс n + 1, указывающий на временной елой, опущен: Qs = (Qn+1)s,
Оставляя только разности первого порядка в конвективной части dRHS/dQ, получим из (13) следующую систему линейных уравнений:
(д7кт + ¿Rt) V + Ci+l/2Ai+l/2 + С+1/2ai_i/2 + Ct+1/2Aj+i/2 + + C+_1/2Aj-1/2 + С —+1/2 Afc+1/2 + C+_i/2Afc _ 1/2 As+1Q =
" -"hs (i4)
где As+1Q = Qs+1 — Qs, С^+1/2 — компоненты расщепленных по знаку собственных значений матриц Якоби потоков через грань m + 1/2,
2.4. Методы решения линеаризованной системы уравнений 2.4.1. Метод приближенной LU-факторизации
В работе [9] для решения линеаризованной системы уравнений применяется метод приближенной LU-факторизации, Для построения факторизованного оператора система (14) переписывается в виде
B + C_+1/2T+ — C+_1/2Ti + С_+1/2Tj+ — C+-1/2Tj +
+ С -+1/2T+ — C+_1/2Tfc_)As+1Q = R(Qs), (15)
где — оператор сдвига па один узел вперед (+) или назад (—) то индекcv m, В = (— RT + ——R*^!7 + С+ /о — С -, 1 /о + Cj, /0 — С •,! /о + С t_i /о ~~ С.
\Ат 2At ) ~1/2 i+l/2 j —1/2 J+l/2 fc—1/2 fc+1/2'
R(Q-) = -R'3Q,-4QA;+Q""+RHS-
и приближенно LU-факторизуетея следующим образом:
B — C+_ 1/2Ti — C+_1/2Tj — C+_1/2Tfc ) B _1 X X (B + C-/2T+ + C _+1/2T/+ — C_+1/2T+) As+1Q = R(Qs). (16)
Система (16) разрешается последовательно в два этапа, на первом из которых вычисления Д5+1О ведутся в направлении возрастания индексов:
Д*+1/2О^ = В-1 [и(0*) + С++— 1/2Д*+1/2О^ +
+ С++— 1/2Д*+1/2Оч—1к + С+— ^Д^1^^— 1], (17)
а на втором — в направлении убывания индексов:
Д-^О^ = Д*+1/2О^ — В—1
Сг+1/2Д5+1^+Ук +
(18)
+ С./+1/2Д-+1Сг7'+1к + Cfc+1/2ДS+lQijfc+1
Значения вспомогательных величин Д-+1/2Сук и приращений Дs+1Qijfc в фиктивных ячейках полагаются равными нулю. Таким образом, в работе [9] используется явная процедура реализации граничных условий,
2.4.2. Метод Гаусса — Зейделя
Для точного разрешения линеаризованной системы уравнений (14) применим итерационный метод Гаусса —Зейделя, заключающийся в том, что решение СЛАУ Д^^ = Др+^ находится в процессе итераций, каждая из которых состоит из двух шагов:
В — С++— 1/2Т— — С+—1/2 — С+— 1/2ТТ) д^1^ + С—+1/2Т+ + С7+1/2Т+ + С —+1/2Т+) д^ = ад-), (19)
В + С—+1/2Т+ + С—+1/2^+ + С—+1/2Т/+) ДP+lQ — — ( С++— 1/2Т— + С+— 1/2Т,— + С+— 1/2Т— ) Д^1^ = И^-). (20)
В качестве начального приближения берется нулевая невязка Д^ = 0, Тогда на первой итерации уравнения (19) совпадают с (17) и из них же и (20) следуют уравнения
В + С—-+1/2Т+ + С-+1/2Т+ + С—+1/2Т+) д^ — ВД^1^ + И^-) = И^-),
из которые, если на них подействовать слева оператором ^В — С+1/2Т" — С+—1/2Т/— — С+—1/2Т— В—1, получается Ш-факторизованная система (16), Таким образом, на пер-
вой итерации метод Гаусса — Зейделя совпадает с методом приближенной Ш-фактори-зации.
Исключая из (19) с помощью (20) величину Д^^^, получим соотношение
В
В + С—+1/2Т+ — С+1/2Т + С—+1/2^+ — С+—1/2^ + С—+1/2Т+ —
— С+—^Т— Д^ — И^-)
1/2
х (С+—1/2^" + С+—1/2Т,— + С+—1/2Т— ) (Д^ — Д^). (21)
С—+1/2Т+ + С—+ 1/2"^?+ + С —+1/2^+ I Х
Видно, что при || Ap+1Q — ApQ|| ^ 0 (21) является решением линеаризованной системы уравнений (15), Таким образом, решение считается найденным, если ||Ap+1Q — ApQ|| < etol, где etoi — заданная точность, или выполнено заданное число итераций iter, достаточное для получения решения с требуемой точностью,
2.4.3. Методы библиотеки Intel MKL
Характерная размерность матрицы СЛАУ (14) для задач моделирования течений в ГТ составляет ~ 105, Матрица является разреженной: число ненулевых элементов составляет ~ 0.001 % от общего количества элементов матрицы. Поэтому для получения точного решения СЛАУ, помимо метода Гаусса —Зейделя, были выбраны методы библиотеки Intel MKL для разреженных матриц: FGMRES и PARDISO,
Для экономии оперативной памяти и уменьшения числа арифметических операций разреженные матрицы СЛАУ принято хранить в специальном формате: матрица представляется в виде одномерного массива, формируемого из ее ненулевых элементов, и двух одномерных массивов, содержащих информацию о координатах расположения ненулевых элементов в исходной матрице. Следует отметить, что в представленном в работе алгоритме от итерации к итерации структура матрицы СЛАУ не меняется, поэтому ее достаточно вычислить один раз — при запуске расчета, а на каждой итерации обновлять только меняющиеся в ходе расчета элементы матрицы. Методы библиотеки MKL работают с форматом хранения разреженных матриц CSE,
FGMRES - метод обобщенных минимальных невязок. Это итерационный метод, основанный на проектировании решения СЛАУ на последовательность подпространств Крылова, Подпространство Крылова представляет собой линейную оболочку векторов из Km(r, A) = spanjr, Ar, A2r,..., Am_ где A — квадратная матрица размерности N X N.
Для решения СЛАУ Ax = b применяется следующий подход. Рассматривается последовательность подпространств Крылова Km(r0, A),m = 1, 2,..., где ro = Axo — b, x0 — начальное приближение. Очередное приближение к решению xm ищется па основе минимизации эвклидовой нормы невязки ||Axm — b|| на текущем подпространстве Крылова Km(r0, A),
В общем виде алгоритм метода FGMRES может быть записан следующим образом:
xo
2) последовательно строится ортонормированный базис в Km (это требует m умножений матрицы па вектор) v1, v2,..., vm;
3) вычисляются коэффициенты а1, а2,..., am так, чтобы вектор xm = x0 + a1v1 + a2v2 + ... + amvm имел минимальную норму.
Шаги 2, 3 повторяются до тех пор, пока ||Axm — b|| > to/||Ax0 — b||.
Очевидно, что вычисление каждого последующего приближения к решению требует большего числа арифметических операций. Для задач, в которых число шагов FGMRES для достижения заданной точности велико, имеет смысл использовать перезапускаемый FGMRES, В этом случае делается заданное число шагов, после чего метод запускается с полученного приближения. Этот процесс повторяется до сходимости с заданной точностью.
Известно, что скорость сходимости итерационных процессов зависит от спектраль-
Ax =
b, а эквивалентную ей MAx = Mb, в которой матрица M в некотором смысле близка к
исходной (АМ—1 близка к единичной) и операция ее обращения реализуется экономично, Такой подход называется предобуелавливанием. Минимизация затрат процессорного времени на расчет не была главной целью данной работы, поэтому для расчетов применялся FGMR.ES без предобуелавливания, В дальнейшем предполагается рассмотрение FGMR.ES с предобуелавливанием,
В методе I 'А Н 1)1Ж) заложена возможность решать СЛАУ как прямым методом (1.Г-разложение), так и итерационным (('С — для симметричных и ('С;Я — для несимметричных матриц), В данной работе использовались оба варианта. Применение итерационного метода актуально для решения задач, в которых необходимо получить набор решений последовательности систем уравнений, незначительно отличающихся друг от друга. Это имеет место при решении нестационарных задач динамики жидкости методом, описанным в разделах 2,1 — 2,3,
При запуске РАШЛЯО с опцией использования итерационного метода для первой решаемой СЛАУ вычисляется ее Ы1-разложение, а все последующие СЛАУ решаются итерационным методом с предобуелавливанием, где в качестве предобуелавливателя используется Ы1-разложение, полученное на первом шаге. Если за 150 итераций решение с заданной пользователем точностью ¿о/ найти не удается, та же система решается прямым методом. Максимальное число итераций фиксировано в РАШЛЯО из тех соображений, что на их выполнение затратится меньше половины времени, требуемого для вычисления Ы1-разложения, Решение следующей СЛАУ ищется вновь итерационным методом с предобуелавливанием, где в качестве предобуелавливателя берется Ы1-разложение, найденное на предыдущем шаге, и т, д,
2.4.4. Реализация краевых условий
Поскольку на (5 + 1)-м слое схема (14) имеет трехточечный шаблон по каждому пространственному направлению, для замыкания системы уравнений (14) при неявном задании краевых условий необходимо ввести дополнительные соотношения только для ближайшего к границе области ряда фиктивных ячеек. Ячейки нумеруются от 1 до Кш1, ш = 1, 2, 3, для направлений г, к соответственно, Кш2 = Кш1 — 1, Кш3 = Кш2 — 1. Рады фиктивных ячеек идут иод номерами 1, 2 и Кш2,Кш1,
Рассмотрим входную границу. На ней фиксирован расход. Это осуществлено следующим образом: поля скоростей в слоях фиктивных ячеек 1^'к и 2^'к определяются так, что расход и циркуляция потока жидкости в этих областях соответствуют заданным. Давление в фиктивные ячейки экстраполируется изнутри области так, чтобы в них сохранялся градиент давления, имеющийся внутри области. Таким образом, на этой границе должно выполняться соотношение
где Ех = сМа§ (1, 0, 0, 0), Р = сМа§ (—(1 + Р23/Р34), 0, 0, 0) Р = сМа§ (р2з/рз4, 0, 0, 0), Р23 и р34 _ расстояния между центрами 2 и 3 и 3 и 4 ячеек соответственно.
На выходной границе фиксировано давление и экстраполируется скорость с первым порядком. На этой границе должно выполняться соотношение
ЕхД-+^ к + РД'+^к + к = 0,
(22)
Е2Д-+1°К12^ + Б^+^к^Тг = °
где Е2 = сИа§ (0,1,1,1), Б = сИа§ (0, —1, —1, —1),
В случае невязких течений на твердой стенке ставится условие равенства нулю нормальной компоненты скорости. Для обеспечения этого условия в фиктивных ячейках сетки на данном типе границ скорости определяются симметрично скоростям из соответствующих ячеек расчетной области относительно границы, а давление — из центрально-разностной аппроксимации соотношения др/дп = 0, Таким образом, на твердой стенке должно выполняться соотношение вида
ЕД-^д*« + мд-^д
¿,3
0,
(24)
где
Е = ^(1,1,1,1), М
2пХ - 1
2П'х Пу
2Пх "Пу
2п2 - 1
2П'у
0
2пт гп 7
\
2П'у >2
-1 0 0 0 0
V 0 2пхП7 2пуп 2п2 — 1 / , п — компоненты нормали к твердой стенке. После расчета поля течения на (5 + 1)-слое в фиктивных ячейках проводится корректировка давления так, чтобы выполнялось соотношение пУр = и(иУ)п+п 1, которое получается из дифференциальных уравнений импульса при условии и • п = 0,
На границе цикличности значения зависимых переменных в фиктивных ячейках получаются с помощью поворота векторов скорости и давления из ячеек расчетной области вокруг оси РК на угол периода а, который определяется количеством лопаток в НА или лопастей в РК, Таким образом, на границе периодичности в случае совпадения оси г с осью РК должно выполняться соотношение типа
ЕД-+1дг,2 + Б^+^^З = 0,
(25)
где
Бс
1 0 0 0
0 сое а — вт а 0
0 вт а сое а 0
0 0 0 1
У расчетных областей имеются границы, отделяющие их от соседних областей-блоков, Фиктивные ячейки каждого из таких блоков совпадают с крайними расчетными ячейками соседнего блока. Данные границы являются искусственными в области течения, Параметры в ячейках с обеих ее сторон связаны только основными уравнениями, и неявная реализация условий на данных границах означает лишь совместное решение объединенной СЛАУ этих блоков, что невозможно в силу чрезмерного увеличения размерности матрицы объединенной СЛАУ, Поэтому на искусственных межблочных границах остается явная процедура заполнения фиктивных ячеек: они заполняются параметрами потока из соседнего блока после расчета внутри него при нулевых невязках в фиктивных ячейках.
2.5. Блочное распараллеливание
Ввиду необходимости при моделировании течений в проточном тракте ГТ разбиения расчетной области на отдельные блоки (см, раздел 1,3) наиболее простым в реализации способом распараллеливания вычислений в данном случае является геометрическое распараллеливание: расчеты каждой итерации метода установления во всех блоках проводятся одновременно на отдельных процессорах, после чего между смежными блоками
происходит обмен граничными значениями и итерации повторяются. Очевидно, что хорошее ускорение такой алгоритм позволит получить, если число ячеек в каждом из блоков одинаково. Данный алгоритм распараллеливания вычислений был реализован для циклической постановки с использованием коммуникационной библиотеки MPI.
3. Результаты расчетов
В данном разделе приводятся результаты расчетов в проточном тракте радиально-осевой гидротурбины РО-910 ГЭС Платановрисси и сравнение их с экспериментальными данными. Оцениваются и сопоставляются скорости сходимости решений, получаемых рассмотренными методами.
3.1. Расчеты стационарного течения в конусе отсасывающей трубы
С целью выяснения влияния выбора метода решения СЛАУ и постановки краевых условий на число итераций, требующихся для сходимости к решению, и на время счета были проведены расчеты невязкого стационарного течения в конусе, подобном конусу ОТ (рис. 2). С методом приближенной LU-факторизации сравнивались следующие подходы: нефакторизованная схема с явными краевыми условиями с применением для решения СЛАУ методов MKL (FGMRES, PARDISO) и метода Гаусса Зейделя; нефакторизованная схема с неявными краевыми условиями с решателем СЛАУ PAR DISO. Поскольку при решении данной задачи число внутренних итераций FGMRES на каждом шаге метода установления не превышало 300, применялся неперезапускаемый метод FGMRES.
Во входном сечении расчетной области задавалось равномерное распределение поля скоростей uin = (cz, 0, 0). Сетка имела размеры 28 х 20 х 20. Шаг то псевдовремени At во всех расчетах равнялся 10 и был выбран из тех соображений, что метод приближен-
-5
-10
lg(div)
500 1000 1500 2000
Рис. 2. Расчетная область
Рис. 3. Истории сходимости к решению: сплошная линия — схема приближенной ЬЦ-факторизации, пунктир — нефакторизованная схема с явными краевыми условиями, штрихпунк-тир — нефакторизованная схема с неявными краевыми условиями
ной Ы1-факторизации при данном значении шага дает решение за минимальное число итераций. Решение вычислялось с точностью &у(и) < 10 -3.
Расчеты показали, что для получения решения нефакторизованной схемы с явными краевыми условиями потребовалось в 1.72 раза, а с неявными краевыми условиями в 30 раз меньше итераций по сравнению с методом приближенной Ы1-факторизации (табл. 1). Истории установления решения задачи различными методами до машинной точности приведены на рис. 3. Отбрасывание в уравнении (12) члена Wm+1/2, отвечающего за третий порядок аппроксимации, приводит к согласованию матрицы оператора и правой части, что позволяет получить ускорение по сравнению с методом приближенной Ы1-факторизации при решении нефакторизованной схемы с явными краевыми условиями в 7.3 раза, с неявными краевыми условиями в 106.6 раза (рис. 4). Однако
Таблица 1. Сравнение различных методов решения СЛАУ по скорости сходимости
Схема Число итераций для Время Ускорение по Замедление
установления до 10_3 счета, с числу итераций по времени
Приближенная LU-факторизация 658 13 1 1
Нефакторизованная FGMRES tol = 10"6 383 398 1.72 30.6
Нефакторизованная
с явными краевыми
условиями PARDISO без итераций 383 3841 1.72 295.46
Нефакторизованная (метод Гаусса—Зейделя) iter = 30 383 92 1.72 7.08
Нефакторизованная
с неявными краевыми
условиями PARDISO tol = 10 -6 22 96 30 7.38
а б в
Рис. 4. Истории сходимости к решению при различных порядках аппроксимации правой части уравнения (12): сплошная линия — третий порядок, штрих — первый порядок; а — схема приближенной IX-факторизации; б — нефакторизованная схема с явными краевыми условиями; в нефакторизованная схема с неявными краевыми условиями
и в этом случае уменьшения времени счета на рассматриваемой задаче получить не удалось.
Увеличение шага по псевдовремени Ат при решении нефакторизованных схем не позволило уменьшить число итераций, требующихся для установления.
3.2. Расчеты нестационарного течения в конусе отсасывающей трубы
3.2.1. Влияние шага по физическому времени
При решении нестационарных задач по методике, изложенной в [9], от величины шага по времени АЬ зависит число итераций по псевдовремени, необходимых для установления решения на каждом временном слое до заданной точности: чем больше шаг АЬ, тем меньше требуется итераций (рис. 5, а). При использовании нефакторизованной схемы с неявными краевыми условиями наблюдается аналогичная картина (рис. 5,6), однако в этом случае число итераций существенно меньше. Поскольку при решении нестационарных задач шаг по времени должен быть достаточно мал, применение нефакторизованной схемы с неявными краевыми условиями по сравнению с исходной методикой позволяет значительно уменьшить общее число итераций.
3.2.2. Моделирование прецессии вихревого жгута в режиме неполной загрузки в конусе отсасывающей трубы
Эффективность работы метода с нефакторизованной схемой и неявной постановкой краевых условий исследовалась на задаче невязкого нестационарного течения в конусе, подобном конусу ОТ. Расчеты проводились на сетке размером 58x40x40. Во входном сечении расчетной области задавались компоненты поля скорости (рис. 6), взятые из расчета течения в НА, РК и ОТ в циклической постановке (см. ниже раздел 3.3). Расчеты соответствовали режиму неполной загрузки работы гидротурбины: расход (ц = 0.595 м3/с, частота вращения РК пц = 73.5 об./мин. Открытие лопаток НА А0 = 51.4%. В выходном сечении конуса задавалось постоянное значение давления. Шаг по времени АЬ = 0.01 с, шаг по псевдовремени Ат = 10. Решение устанавливалось до выполнения условия Шу(и) < 10-3.
а б
Рис. 5. Характерные истории сходимости к решению на временном слое: сплошная линия — № = 0.001с, штрих — № = 0.01 с, штрнхпунктир — № = 0.1с; а — схема приближенной Ьи-факторизации, б— нефакторизованная схема с неявными краевыми условиями
Для радиально-осевых ГТ в режиме неполной загрузки при A0 < 65 % характерно наличие прецессирующего вихревого жгута за РК. На рис. 7 показан вихревой жгут, а — полученный в эксперименте, его визуализация была осуществлена путем вдувания воздуха в поток жидкости и зарисовки места скопления пузырьков, б — расчетный, визуализированный поверхностью постоянного давления.
Пульсации давления в точке 1 на стенке конуса (рис. 7, а), обусловленные наличием прецессирующего жгута, а также их спектральный анализ приведены на рис. 8. Отношение частоты прецессии вихревого жгута / к частоте вращения РК fR, найденное расчетным путем, составило 0.2167, в эксперименте для этого отношения получено значение 0.23.
При расчете 1000 шагов по времени (10 с) решение нефакторизованной схемы с неявными краевыми условиями потребовало в 55 раз меньше итераций метода установления, чем решение приближенно LU-факторизованной схемы. Однако для этого при использовании для решения СЛАУ итерационного метода PARDISO (tol = 10-6) потребовалось в 3 раза больше машинного времени. Следует отметить, что для рассматриваемой задачи характерно существенное замедление скорости сходимости при решении приближенно LU-факторизованной схемы после достижения невязки величи-
Рис. 6. Профили радиальнои С г-, осевой Сг и окружной Си компонент скорости во входном сечении расчетной области
Рис. 7. Вихревой жгут за рабочим колесом: а — эксперимент, б— расчет
а б
Рис.. 8. Пульсации давления в точке 1 (см. рис.7, а) (а), спектральный анализ пульсаций (б): сплошная линия — расчет, штрих — эксперимент
Рис. 9. Характерные истории сходимости к решению на временном слое: сплошная линия — приближенно Ьи-факторизованная схема, штрих — нефакторизованная схема с неявными краевыми условиями
ны 10-4 (рис. 9), в то время как при решениии нефакторизованной схемы с неявными краевыми условиями такой тенденции не наблюдается. Поэтому при расчетах с большей точностью установления использование не4^акторизованной схемы с неявными краевыми условиями позволит уменьшить число итераций в сотни раз.
3.3. Расчеты нестационарного течения в направляющем аппарате, рабочем колесе и отсасывающей трубе. Ускорение на нескольких процессорах
Еще один способ уменьшения вычислительных затрат на проведение нестационарных расчетов в проточном тракте гидротурбины упрощение постановки задачи и применение технологий параллельных вычислений.
Наиболее простой является циклическая нестационарная постановка, которая хотя и не позволяет учитывать окружную неравномерность потока и связанные с ней нестационарные эффекты, как будет показано ниже, дает возможность адекватно модели-
ровать такое важное для проектировщиков явление, возникающее при неоптимальных режимах работы Г Г. как низкочастотные пульсации давления в проточном тракте, вызванные прецессией вихревого жгута.
Еще одно достоинство циклической нестационарной постановки состоит в том, что шаг по времени At можно задавать довольно большим, что нельзя сказать о полной нестационарной постановке, где для получения правильного описания взаимодействия лопаток НА с лопастями РК он должен составлять величину порядка 1/100 периода вращения РК.
Были проведены расчеты невязкого и вязкого нестационарных течений в проточном тракте ГТ в режиме неполной загрузки Qn = 0.595 м3/е при пц = 73.5 об./мин, A0 = 51.4% в циклической и полной постановках. В первом случае расчетная область состояла из пяти блоков: одного межлопаточного канала НА, одного межлопастного канала РК, одного конуса и двух выходных диффузоров ОТ (см. рис. 1,6). Сетки имели размеры 46 х 17 х 19 46 х 17 х 15 26 х 40 х 40 и 40 х 25 х 15 соответственно. В полной постановке расчетная область состояла из 23 блоков: 16 межлопастных каналов РК, четырех секторов конуса ОТ, колена ОТ и двух выходных диффузоров ОТ (см. рис. 1, о). Сетки имели размеры 48 х 17 х 15 26 х 24 х 24, 9 х 40 х 40 и 35 х 25 х 15 соответственно. At
Для моделирования вязкого течения была выбрана k — е-модель турбулентности Кима Чена, поскольку при использовании стандартной k — е-модели в диффузоре ОТ возникли зоны рециркуляционного течения (рис. 10). На рис. 11 приведены мгновенные структуры течения за РК, полученные в расчетах невязких и вязких течений. В расчетах, как и в эксперименте, наблюдалась ярко выраженная прецессия вихревого жгута. Расчеты показали прецессию вихревого ядра с периодами 3.6-4 с. Отношение частоты прецессии вихря f к частоте вращения РК fR составило 0.2-0.22 (табл. 2), по экспериментальным данным эта величина равна 0.23. Сравнение рассчитанных спектров пульсаций давления в точке 1 (см. рис. 7, а) с полученным в эксперименте приведено на рис. 12. Видно, что по форме вихря и частоте прецессии расчеты показали хорошее согласование с экспериментом.
Расчеты в циклической постановке проводились на пяти процессорах, в полной постановке — на шести процессорах на кластере информационно-вычислительного центра НГУ, в состав которого входили процессоры Intel Xeon 2,6GHz Quad-Core. В полной постановке течение в РК рассчитывалось на четырех процессорах, в конусе ОТ — на одном, в колене и выходных диффузорах ОТ — также на одном. В полной постановке распределение числа обрабатываемых ячеек по процессорам было следующим: 4 х 48960 — 59904 — 42525, в циклической - 14858 — 11730 — 41600 — 15000 — 15000. В связи с таким распределением время счета одного временного слоя в вязкой циклической постановке составило 1300 с, в полной — 1700 с.
Таблица 2. Отношение рассчитанных частот прецессии вихря к частоте вращения РК
Течение / Постановка Отношение f/fn
Невязкое/Циклическая 0.2027
Невязкое/Полная 0.2180
Вязкое / Циклическая 0.2190
Вязкое / Полная 0.2170
Эксперимент 0.2300
Рис. 10. Меридиональное сечение диффузора; расчет по стандартной к — е-модели турбулентности
I
II
а б
Рис. 11. Вихревой жгут. Расчеты невязкхнх) (I) и турбулентшнх) (II) течения; а циклическая, б полная постановка
ю1 //Л
п
Ю1 ///*
Ю1 //л
ю1 ///я
Рис. 12. Спектральный анализ; расчеты невязксих) (I) и турбулентнсих) (II) течения; сплошная линия расчет, пунктир эксперимент; а циклическая, б полная постановка
В силу того что число ячеек в конусе ОТ существенно превышало число ячеек в остальных блоках, ускорение па пяти процессорах в циклической постановке составило всего 2.2 раза. Ускорение па пяти процессорах расчета в циклической постановке в данном случае составило 2.36 раза.
Заключение
Дня улучшения сходимости итераций в нестационарных задачах построена схема с неявной реализацией граничных условий. В отличие от базового метода решения нестационарных трехмерных уравнений движения несжимаемой жидкости |9| решение СЛАУ, возникающей при неявной аппроксимации уравнений движения жидкости с внесенными в нее соотношениями па границах области, находилось точно с помощью итерационного алгоритма PARDISO из стандартной библиотеки Intel Fortran MKL. Получено ускорение сходимости итераций в сотни раз, причем с повышением точности установления
эффект ускорения увеличивался. Однако точное обращение матрицы СЛАУ по сравнению с факторизованной требует в 80 раз больше процессорного времени. Поэтому для получения существенного уменьшения временных затрат на проведение расчетов в перспективе необходимо сократить время обращения матрицы.
Список литературы
fl] Harlow F., Welch J. Numerical calculation of time-dependent viscous incompressible flow with free surface // Phvs. Fluids. 1966. Vol. 9, No 12. P. 2182-2189.
[21 Amsden A., Harlow F. The SMAC Method. Los Alamos Sci. Lab. Rep. No LA-4370. Los Alamos, 1970.
[3] Белоцерковский О.M. Численное моделирование в механике сплошных сред. М.: Наука, 1984.
[4] Рыков В.В. Численное моделирование нестационарных течений несжимаемой жидкости // Журн. вычисл. математики и матем. физики. 1985. Т. 25, № 5. С. 789-793.
[5] Владимирова Н.Н., Кузнецов Б.Г., Яненко Н.Н. Численный расчет симметричного обтекания пластинки плоским потоком вязкой несжимаемой жидкости // Некоторые вопросы прикладной и вычислительной математики. Новосибирск: ВЦ СО АН СССР, 1966. С. 186-192.
[6] Chorin A.J. A numerical method for solving incompressible viscous flow problems //J. Сотр. Phvs. 1967. Vol. 2. P. 12-26.
[7] Chakravarthy S.R., Osher S. A new class of high resolution TVD schemes for hyperbolic conservations laws // AIAA Pap. 85-0363. 1985.
[8] Rizzi A.W., Eriksson L.E. Computation of inviscid incompressible flow with rotation // J. Fluid. Mech. 1985. Vol. 153, No 12. P. 275-312.
[9] Чёрный С.Г., Чирков Д.В., Лапин В.Н. и др. Численное моделирование течений в турбомашинах. Новосибирск: Наука, 2006.
[10] Chen Y.S., Kim S.W. Computation of Turbulent Flows Using an k-e Turbulence Closure Model. NASA CR-179204.
Поступила в редакцию 28 февраля 2011 г.
