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

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

Том 17, № 6, 2012

Численный алгоритм моделирования

пространственных течений несжимаемой жидкости

>к

на подвижных сетках*

А. Ю. АвдюшЕнко, С. Г. Чёрный, Д. В. Чирков Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected], [email protected]

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

Ключевые слова: нестационарные течения, движущиеся сетки, условие геометрической консервативности, гидротурбина.

Введение

В [1] для моделирования установившихся (стационарных или периодически нестационарных) течений в проточном тракте (ПТ) гидротурбины (ГТ) при неподвижных лопатках направляющего аппарата (НА) предложен численный метод решения трёхмерных уравнений несжимаемой жидкости, хорошо зарекомендовавший себя в работе на фиксированных сетках. В случае переходных течений, возникающих вследствие изменения положения лопаток НА, требуется обобщение метода на подвижные сетки, подстраивающиеся со временем под положения жёстких границ расчётных областей. Выписать основные уравнения движения несжимаемой жидкости в криволинейных координатах, зависящих от времени, не представляет трудности. Однако проведение правильной дискретизации этих уравнений для построения эффективного численного метода требует выполнения определённых условий. Главным из них является условие геометрической консервативности (УГК), которое заключается в сохранении без какого-либо возмущения постоянного однородного потока при расчёте его на двигающейся сетке.

В настоящей работе метод [1] обобщается на задачи с подвижными сетками. Параллельно с проводимой дискретизацией уравнений рассматриваются известные подходы построения численных методов на подвижных сетках, предложенные в [2-11]. Осуществляется их сравнительный анализ и сопоставление с предлагаемым авторами методом.

* Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00475-а) и Интеграционного проекта СО РАН № 130.

1. Основные уравнения

1.1. Дифференциальная форма записи в декартовой системе координат

Течение вязкой несжимаемой жидкости в проточном тракте ГТ, как ив [1], описывается нестационарными трёхмерными уравнениями Рейнольдса, имеющими в декартовой системе координат (х, у, г) = (х\, х2, х3) следующий вид:

Зщ дх,

0,

дт. ~д1

г дтгт, др д

дх

,

д хг д х,

Vей

дт + дт'

дх, дх

+ Iг

2

1, 2, 3,

где т1, т2, т3 — компоненты вектора скорости; р = рс + -к; рс — гидростатическое

3

давление, делённое на плотность жидкости; к — кинетическая энергия турбулентности; I = х^2 + 2т2ш; I = х2ш2 — 2т1ш; 13 = д, д — ускорение свободного падения; ш — угловая скорость вращения рабочего колеса вокруг оси х3; ш = 0 в остальных элементах проточной части.

Величина veff есть сумма молекулярной V и турбулентной V* вязкостей

Veff = V + VI.

(3)

Турбулентная вязкость V* определяется по стандартной к — е-модели турбулентности.

Для удобства дальнейшей работы с уравнениями запишем их в виде векторного уравнения

+ Ех + Оу + Н = Г,

где д = (р,т1,т2,тз)Т, К* = diag (0,1,1,1), Г = (0, Ь, 12, 1з)

т

Е

( т1 \

т2 + р — тЦ

^1^2 — Т21

\ ^1^3 — Т31 /

О

( ^2 \

^1^2 — Т12

Щ + р — Т22 \ ^2^3 — Т32 / дтг дт,

Н

( т \

^1^3 — Т13 ^2^3 — Т23

\ Щ + р — Т33 )

Т

г,

Veff

+

дх, дх.

(4)

(6)

Перейдём от дифференциального векторного уравнения (4) к интегральному, выполняющемуся на произвольном подвижном объёме.

1.2. Интегральные законы сохранения для подвижного объёма

Интегрируя уравнение (4) по произвольному подвижному объёму V(¿) с границей дУ (¿), ориентация которой задаётся внешней единичной нормалью п, и применяя теорему Гаусса — Остроградского, получим

К* У дд¿V + У КЖ = I ГёУ, (7)

V (*) ЭУ(Ь) У(Ь)

где = п^б"; dS - элемент поверхности дУ(¿); матрица К имеет структуру К = (Е, С, Н). Интеграл в первом слагаемом левой части уравнения (7) преобразуется по правилу Лейбница для производной от интеграла с переменными пределами [12]

9 Q ™ 9

ltdV = dt

d x

QdV — ф Q— ■ dS,

V(t) V(t) dV(t)

где x — радиус-вектор элемента поверхности dS;

d x

—- ■ dS — Xt ■ dS — (xt, yt, Zt) I dSy

0t 1 dS.

dSx dSy I —

нормальная скорость перемещения этого элемента. Подставляя (8) в (7), получим

Б/д У QdУ + £ (^Я - Б^х, ■ dS) = I FdУ.

v (t)

dV (t)

V(t)

Преобразуем второй член в круглых скобках (9) под поверхностным интегралом

RtQxt ■ dS — Rt

Rt

( pxt ■ dS \

w\xt ■ dS w2xt ■ dS \ w^xt ■ dS J ( pxt pyt pzt \ WiXt Wi yt WiZt W2Xt W2 yt W2Zt \ W3Xt W3yt W3Zt J

Rt

dSx dSy

dSz

( p (XtdSx + ytdSy + ZtdSz) \ Wi (xtdSx + ytdSy + ZtdSz) W2 (xtdSx + ytdSy + ZtdSz) У W3 (XtdSx + ytdSy + ZddSz) J / 0 0 0 \

WiXt Wiyt WiZt W2Xt W2yt W2Zt \ W3Xt W3yt W3Zt J

dSx dSy

dSz

сгруппируем его с матрицей потоков К и обозначим объединённую матрицу через К*

Kt

Wi

W2 + p — Til — Wi Xt W1W2 — T2i — W2Xt

\ WiW3 — T3i — W3Xt

W3

W1 W3 — T13

W2W3 — T23 2

\

W2

W1 W2 — T12 — W1 yt W2 + p — T22 — W2yt W2W3 — T32 — W3yt w2 + p — T33 — W3Zt J

W1 Zt W2Zt

:iq)

Окончательно интегральные законы сохранения массы и количества движения для подвижного объёма У (¿) примут вид

Rtd i QdV + ф KtdS

FdV.

V(t)

dV (t)

V(t)

1.3. Условие геометрической консервативности в интегральной форме

Уравнения, описывающие течение жидкости, должны сохранять со временем постоянный однородный поток в отсутствие внешних сил. Поэтому равенства F — 0, Q — const

и условие замкнутости у) вБ = 0 выделенного объёма V(¿), подставленные в уравне-

дУ (*)

ние (11), дают тождественное скалярное уравнение

ду(I) = / х* ■ ¿Б, (12)

дУ (*)

которое также является следствием правила Лейбница (8) и называется условием геометрической консервативности. Оно имеет простую геометрическую интерпретацию — изменение выделенного объёма во времени равно сумме объёмов, образованных геометрическими элементами поверхности объёма при их движении. Впервые понятие УГК было введено в [2] для конечно-разностного метода, а позднее распространено на метод конечных объёмов [3, 4]. В работе [4] было предложено вычислять скорости движения граней ячеек исходя из известных перемещений узлов так, чтобы условие геометрической консервативности автоматически выполнялось. При этом в [4] рассмотрен только случай двумерных течений. В настоящей работе указанный подход получения метода, удовлетворяющего УГК, распространён на случай трёхмерных течений.

В [4-6] показано, что нарушение УГК приводит к возникновению осцилляций решения, вызываемых только движением сетки. В работах [6, 7] предложены удовлетворяющие УГК схемы соответственно первого и второго порядка по времени для двумерного и трёхмерного случаев движения сетки.

Если в дискретном аналоге уравнения (11) не удовлетворить условие геометрической консервативности, то на задаче с решением в виде постоянного однородного потока численный метод на подвижной сетке даст возмущённое, отличное от начальных данных решение. Ниже в разделе 2 проводится обобщение метода [1] на случай подвижных сеток. Для полноты изложения приведём также геометрический закон сохранения в дифференциальной форме.

1.4. Условие геометрической консервативности в дифференциальной форме

Наряду с физической областью (х,у,г) с подвижной сеткой и подвижным элементом объёма физического пространства вУ рассмотрим вычислительную область (£,г/,() с неподвижным элементом объёма вычислительного пространства вУ0. Между декартовыми координатами физической области и криволинейными координатами вычислительной области устанавливается взаимно-однозначное отображение

£ = C(x,У,z,t),

П = п(х, У, г, г),

С = С (x,y,z,t),

г' = г. (13)

Рассмотрим матрицу

А

/ £х £у \

Пх 'Пу Vz V*

Сх Су Су С*

V 0 0 0 1 /

14)

и определим величину

3 = ёе! А. (15)

Тогда связь между элементами объёмов в физическом и вычислительном пространствах задаётся равенством

3 ■ ¿У = ¿Уо. (16)

Следует отметить, что сетка и объёмы в вычислительном пространстве фиксированы и не зависят от времени. Но из-за подвижности сетки и объёмов в физическом пространстве преобразование координат (13) и соответствующая ему матрица Якоби (14) зависят от времени. В силу теоремы Гаусса — Остроградского УГК (12) переписывается в виде

¿V = у а1у (х*) ¿V, (17)

а с учетом (16) и независимости ¿У0 от времени —

/ д (1) ¿У° = / 3^ (х*) ¿Уо. (18)

Vo Vo

Из (18) следует дифференциальная форма геометрического закона сохранения

д 1

З-3 - х* = 0. (19)

2. Численный метод

2.1. Дискретизация уравнений для подвижного объёма

Расчётная физическая область разбивается на элементарные ячейки Упк в виде криволинейных шестигранников. Центрам ячеек приписываются осреднённые по объёмам значения переменных и массовых сил

(ЪчкУг3к)П = У ^¿У, (¥гзкУг]к)П = У Fn¿V. (20)

V" V.",

гук гук

Неявная конечно-объёмная аппроксимация уравнения (11) на ячейке У^к дает

з(дУГ1 - 4(дУ)п + (дУ)п-1 _ ,

-' 1 ^ ^— = яш* , (21)

2Аг V ; ' V у

где ДЬ — шаг по времени, п — номер слоя по времени, Уп — объём ячейки на п-м слое по времени. Правая часть (21) имеет структуру

ЯШ* = - ((К*8)г+1/2 - (К*Я)^-1/2 + (К^)^ - (К*8)^-1/2 + + (К*Я)к+1/2 - (К*Я)к-1/2) + ~РУ,

где выражения (К*Я)^+1/2 , (К*Я)^-+1/2 , (К*Я)к+1/2 представляют собой разностные потоки через грани Si+1/2, Sj+1/2, Sk+1/2 ячейки с номером %]к и объёмом У^к. Определим векторы

^+1/2 = ^, ^+1/2 , ^'+1/2 = (Sx, Sу, )j+l/2 ,

Як+1/2 = (5Х, Sу, )к+1/2 (22)

как нормали к граням Si+1/2, Sj+1/2, Sk+1/2 ячейки с номером %г]к. Длины этих векторов равны площадям соответствующих граней данной ячейки. Направления нормалей выбираются таким образом, чтобы к граням Si+1/2, Sj+1/2, Sk+1/2 они были внешними, а к граням Si-1/2, Sj-1/2, Sk-1/2 — внутренними по отношению к ячейке сетки с индексом %г]к. Такой выбор нормалей упрощает их построение, поскольку при этом выделяются одни глобальные для всех ячеек сетки, положительные по каждой из координат направления. Нормали к граням ячейки определяются таким образом, чтобы выполнялось соотношение

Si+i/2 — Si-i/2 + Sj+i/2 — Sj-i/2 + Sfc+i/2 — Sfc-i/2 — 0, (23)

являющееся следствием интегрального условия замкнутости j) dS — 0.

dv (t)

2.2. Дискретное условие геометрической консервативности

Подставив в (21) значения F — 0, Q — const и учитывая дискретное условие замкнутости ячейки (23), получим дискретное условие геометрической консервативности

3V n+i — 4Vn + Vn-i

— — (xt ■ S)i+i/2 — (xt ■ S)i-i/2 + (xt ■ S)j+i/2 — (xt ■ S)j-i/2 +

+ (xt ■ S)fc+i/2 — (xt ■ S)fc-i/2 . (24)

Соотношение (24) связывает способы вычисления объёма ячейки и нормальной составляющей скорости движения грани xt ■ S.

2.3. Скорость движения грани ячейки

Выведем из (24) конкретные выражения для скорости движения грани ячейки xt ■ S.

2.3.1. Одномерный случай

В одномерном случае скорости движения граней и нормали вырождаются в скалярные величины и УГК (24) принимает вид

3Уп+1 - 4УП + уп-1

—-2А; + i = (*^+1/2 - (*Ж_1/2 . (25)

Объёмы ячеек вычисляются как У¿ = £¿+1/2 — 1/2, где £¿±1/2 — координаты краев 1-й ячейки. Левая часть уравнения (25) может быть преобразована как

ЗУТ+1 - + V"-1

2ДЬ

3 (£¿+1/2 — £¿—1/2) "+1 — 4 (£¿+1/2 — £¿-1/2)" + (£¿+1/2 — £-1/2)'

2ДЬ

3гГп+1 _ 4£п + xn— 1 о£п+1 _ 4£п I ^га-1

3£+1/2 4£+1/2 + £+1/2 3£—1/2 4£—1/2 + 1/2

2ДЬ 2ДЬ

/ Зжга+1 — 4£п + £п—1 \ / Зжга+1 — 4£п + £п—1

V 2Д / ¿+1/2 V 2Д /¿—1/2

(26)

Из (26) следует, что для тождественного выполнения УГК (25) скорости движения граней ячеек в одномерном случае должны рассчитываться по направленной назад разности второго порядка

Зж^1 - 4£п + £п—1

£ =-х-. (27)

4 2Дг v 7

2.3.2. Двумерный случай

В двумерном случае объёмы ячеек есть их площади, грани ячеек — отрезки прямых. За основу построения выражения для скорости движения грани ячейки в этом случае взята методика, предложенная для двуслойной схемы в работе [4]. Для двуслойной схемы УГК имеет вид

уп+1 _ уп

¿7 Д ¿7 = (X ■ й^+1/2 — (х ■ S)¿—1/2 + (х ■ 8)^+1/2 — (х* ■ 8)^—1/2, (28)

что можно переписать как

ДУП

¿77

"д^ = (х* ■ ^¿+1/2 — (х* ■ 1/2 + (х* ■ ^+1/2 — (х* ■ 8)^ —1/2 , (29)

где

ДУп = У¿n+1 — У¿n (30)

есть изменение площади ячейки при переходе с п-го слоя на п + 1-й слой. В работе [6] показано, что

ДУП = У+1/2 — уг—1/2 + У,п+1/2 — уп—1/2, (31)

где

Ут±1/2 = (хп+1 — хп) 8 2+ 8 , т = г,],

V 2 /т±1/2

представляют собой площади фигур, образованных гранями-отрезками £т±1/2 при их движении за время ДЬ. Плюсы при У+1/2, У+1/2 и минусы при У—1/2, У-1/2 в правой части (31) обусловлены тем, что направления нормалей к отрезкам £¿+1/2, £7+1/2 являются внешними, а к £¿^/2, £7-1/2 — внутренними по отношению к ячейке сетки с индексом г].

Подставляя (31) в (29), получим

т/п т/га т/га т/га ^+1/2 ^—1/2 ^+1/2 Vj—1/2 _ - - - —I— - - - -

Аг Аг Аг Аг = (х* ■ й^+1/2 - (х* ■ S)i-1/2 + (х* ■ S)j+1/2 - (х* ■ S)j-1/2 . (32)

Следовательно, если положить

У п

(х* ■ 8)т±1/2 = тА±1/2 , т = г,J, (33)

то УГК будет выполняться точно.

В нашем случае трёхслойной схемы УГК имеет вид

3Уп+1 - 4Уп + Уп—1

—-^-— = (х* ■ - (х* ■ 8)<_1/2 + (х* ■ ^+1/2 - (х* ■ ^—1/2 , (34)

что можно переписать как 3АУп — АУп—1

2А = (х* ■ йХ+1/2 - (х* ■ 1/2 + (х* ■ ^^'+1/2 - (х* ■ 8),7 —1/2 , (35)

где величина АУп определена в (30) и вычисляется также по формуле (31). Подставляя (31) в (35), получим

3Уп — Уп—1 \ / 3Уп _ Уп—/ 3Уп — Уп— П / 3Уп — Уп—1

+

— I - I I - I —

2Аг Л+1/2 V 2Аг Л—1/2 V 2Аг А+1/2 V 2Аг / j—1/2 = (х* ■ S)i+1/2 - (х* ■ S)i—1/2 + (х* ■ ^'+1/2 - (х* ■ —1/2 . (36)

Следовательно, если положить

_ (3УП - УП—1)т±1/2

то УГК будет выполняться точно.

(х* ■ 8)т±1/2 =-2А-, т = г,j, (37)

2.3.3. Трёхмерный случай

В трёхмерном случае авторами развита идея работы [4] для двумерного случая. Дискретное условие геометрической консервативности (24) перепишем как 3А У п А У п— 1

— = (х* ■ йХ+1/2 - (х* ■ 1/2 + (х* ■ й^+1/2 - (х* ■ 8),7 —1/2 +

+ (х* ■ й)к+1/2 - (х* ■ й)к—1/2 , (38)

где А Уп = Уп+1 - Уп. По аналогии с двумерным случаем изменение объёма ячейки представим в виде

п п п п п п п

АУ = *г+1 /2 - 4—1/2 + У^'+1/2 - Vj—1/2 + Ук+1 /2 - Ук—1/2,

где Ут±1/2 — объёмы фигур, образованных гранями $т±1/2 при их движении за время А*. Подставим соотношение (39) в уравнение (38)

от™ _ т/п—1 ЧЛ/П _ Т/п—1 _ Т/п—1 °Л/П _ Т/п—1

3 4+1/2 4+1/2 3 4—1/2 4—1/2 3 ^+1/2 ;/+1/2 3 ^ —1/2 V/ —1/2

+---—Г"---—Г"--Г

2А* 2А* 2А* 2А*

п п— 1 п п— 1

. 3Кк+1/2 4+1/2 3Кк—1/2 Ук—1/2 / Сч / Сч

+-2А*---2А- = (Х ■ Й)г+1/2 - (Х ■ Й)г—1/2 +

+ (Х ■ 8/+1/2 - (Х ■ —1/2 + (Х ■ Я)к+1/2 - (Х ■ Я)к—1/2 . (40)

Из (40) заключаем, что для точного выполнения УГК достаточно положить

(3Уп - Vп—1)т±1/2

(Х ■ й)т±1/2 = -2А*-' т = % (41)

2.4. Метод вычисления объёмов Vп и

обусловливающий точное выполнение равенства (39)

Комплексы (х4 ■ Я)т+1/2, вычисленные по формуле (41), будут давать точное выполнение дискретного УГК (24), если соотношение (39) выполнено точно, для чего необходимо определённым образом вычислять объёмы Ут±1/2, заметаемые гранями ячеек при движении сетки. Очевидно, способ их вычисления зависит от способа расчёта самих объёмов Vп. В [11] сформулированы общие принципы выполнения (39) на основе метода неопределённых коэффициентов, однако предложенная там методика применима только к одномерным и двумерным задачам.

В настоящей работе предлагается следующий метод вычисления объёмов фигур, входящих в равенство (39).

Пусть объём Vink образован вершинами с радиус-векторами

хп хп хп хп

хг,/,к ' хг,/+1,к ' Хг,/+1,к+1 , Хг,/,к+1 ,

хп хп хп хп (42)

Хг+1,/,к ' Хг+1,/+1,к ' Хг+1,/+1,к+1, Хг+1,/,к+1.

Объём ячейки Vijk вычислим как сумму объёмов шести треугольных пирамид, имеющих общее ребро, соединяющее вершины х*/к и хп+1 /+1 к+1 (рис. 1, а):

У/к = 6 (Хп+1 ,/+1, к+1 — Хг,/, к) ( [(хг,/+1, к — Хг,/ , к+1) Х (хг,/+1 ,к+1 — Хг,/, к)] + + [(Хп/,к+1 — Хп+1,/,к) Х (Хп+1,/,к+1 — Хп/,к^ +

+ [(хг+1,/,к — Хп,/+1,к) Х (хп+1,/+1,к — Хг,/,к. (43)

Согласно введённым обозначениям, объём V—1/2, заметаемый гранью 1/2 при ее движении за время А* = ¿п+1 — *п, образован вершинами

хп хп хп хп хп+1 хп+1 хп+1 хп+1 (44)

хг,/,к ' хг/+1,к ' Хг,/+1,к+1 ' Хг,/,к+1 ' хг,/,к ' хг/+1,к ' Хг,/+1,к+1 ' Хг,/,к+1. V44/

Заметим, что обозначения вершин объёма V—1/2 совпадают с обозначениями (42) вершин объёма Угпк с точностью до замены индексов (■)п+1 на индексы (■)п+1 (рис. 1, б).

Рис. 1. Объём "У^к (а) и объём "^—1/2, заметаемый гранью Б—1/2 при движении сетки за время

АЬ = Г+1 - Ьп, (б)

Объёмы У—1/2 нына (■)п+1:

вычисляются по

той же формуле (43), в которой индексы (-)п+1 замене-

У п

Уг—1/2

1 (хп+1 _ х

6 1хад + 1,к+1 Xi

г,з,к) (К^+м- хЪ,к+1) х (хп+1,к+1- хЪ,к)] +

+ [(х1зк+1- ) х (хп+к+1- х.^ +

+ [(хШ - хЬ+1,к) х (х3+1,к - хХзк^ ) •

(45)

Аналогично с использованием формулы (43) вычисляются объёмы Ур—1/2 и Для

нахождения Ур—1/2 индексы (■)п+1 заменяются на (■ )п+1, для нахождения Укг—1/2 — (■)п+1 заменяются на (-)п+1.

С помощью элементарных преобразований показано, что при вычислении объёмов

Уф и У£+ по формуле (43), а объёмов У^±1/2, ^п±1/2, Ук±1/2 п0 формулам вида (45) равенство (39) выполнено тождественно. При проведении преобразований в правых частях (43) и (45) предварительно приведены подобные слагаемые. Так, (43) эквивалентно равенству

У1пк

1

-|х

6

— ь 11) х х

п

■^З^Ъ

(хп+1,7 + 1,к+1 хЬ,к) ( [(хГ!п,з + 1,

+ [(хп.?,к+1 - х+1 ¿,к) х х+1Л,к+1] + п — хп ) у хп

п

ij+1,k+1

] +

+ [(х?+, ,. - хг.

Аналогично упрощается выражение (45).

2.5. Обобщение метода искусственной сжимаемости [1] на подвижные сетки

2.5.1. Метод искусственной сжимаемости

Метод искусственной сжимаемости заключается во введении в уравнение неразрывности производной по псевдовремени от давления, а в уравнения количества движения — производных по псевдовремени от соответствующих компонент скорости. Модифицированное уравнение (11) запишется в виде

птд , Ид

И дТ + Ит

у QdV +

У(4) дУ(4)

К dS

(46)

V (г)

где ИТ"

Кв

diag(1 Д Д Д), ( в^ц

\

т"2 + р — Гц — -Ш1Хг ^1^2 — Т21 — т2Хг ^1^3 — тз1 — т3Хг

вт2

^1^2 — Т12 — т1 уг

т + Р — г22 — т2Уг

т2тз — тз2 — тз уг

вт

^1^3 — Т13 — т1 г

т2тз — Т2з — т2 г т3 + р — тзз — тзгг )

(47)

В (46) в — коэффициент искусственной сжимаемости. Уравнение (46) решается численно. При этом на каждом шаге по физическому времени Ь проводится установление решения по псевдовремени т.

2.5.2. Неявная конечно-объёмная аппроксимация

Дискретизация уравнения (46) даёт

0^п+1, 8+1 _ Qn+1, 5

Ит

Ат

+И

4 3(QV)п+1, 8+1 — 4 (QV)п + (QV)п—1

2АЬ

(ЯШ4)

Пп+1, 8 + 1

(48)

где Ат' АЬ — шаги по псевдовремени и физическому времени соответственно, в — номер итерации по псевдовремени, п — номер слоя по времени. Дискретизация (48) имеет первый порядок аппроксимации по псевдовремени т и второй порядок аппроксимации по физическому времени Ь. Правая часть имеет структуру

ЯШ 4 = —((Кв ■ Я)г+1/2 — (Кв ■ 1/2 +

+ (Кв ■ Б),-+1/2 — (Кв ■ Я)/ —1/2 + (Кв ■ Я)к+1/2 — (Кв ■ Я)к—1/2) + FV' где разностные потоки через грани ячеек представляются в виде суммы конвективного (невязкого) и вязкого потоков:

Кв ^ т+1/2

Кв,гпи ■ т+1/2 + (Кв,тз ■ Я

т+1/2 '

(49)

Из (47) следует, что движение узлов сетки влияет на вычисление только невязких потоков. Невязкие потоки имеют вид

К4 • Я

в, гп^ ^

( ви \

(и — ид) т + pSx (и — ид) т2 + pSy \ (и — ид) тз + pSz У

где и = w • Я, иа = х* • Я. Для вычисления разностных невязких потоков используется МиБСЬ-схема [13], обеспечивающая третий порядок аппроксимации по пространству

К

Б

1

в,т'о ~Jm+ а 2

О [(Кв,т' Ш + Кв, гп'о • Бт+1/2

- |Л| - д

(51)

где

0,Ь 0,т + ^ = Ът+1 — 4

2 4

3 (дт Ът-1) + 3 (Ът+1 Ът)

24

3 (Ът+1 — Ът) + 3 (Ът+2 — Ът+1)

(52)

Под модулем матрицы |Л*| в (51) подразумевается разность

|Л*| = А+ - А-

матриц Л+ и Л- из разложения матрицы Якоби невязкого потока

А(Ъ)

(53)

д(Кв,_ (д) • Б) 0 в$у

/ш1Бх + и - и я

д Ъ + и - и

_ $ у 3 3 + и

на сумму

Л1 = Л+ + Л

(54)

(55)

Расщепление матрицы Л* на сумму матриц Л+ и Л (55) осуществляется на основе наличия у Л* действительных собственных значений

Аи = и - и, Л

3,4

и - 2 и ± с,

(56)

где с

\

и - ) + в ($2 + $2 + $2). Матрицы Л+ и Л имеют неотрицательные

и неположительные собственные значения. Расщепление (55) может быть осуществлено различными способами, рассмотренными в [1]. В этой же работе предложен экономичный метод решения уравнений (48).

ь

я

3. Краевые условия на подвижной твёрдой границе

Для скорости на подвижной границе Г задаётся условие прилипания

w|г = Уг, (57)

где Уг — скорость движения твёрдой границы.

Для выяснения вопроса о необходимости модификации условий для давления на твёрдой стенке в случае её движения рассмотрим простейшую модельную ситуацию. Так как для расчёта давления на поверхности тела используется уравнение количества движения в проекции на нормаль п, то рассмотрим локально одномерную задачу,

в которой твёрдой стенкой является плоскость x = const. Нормальная составляющая уравнения количества движения к этой стенке есть

1 . . ut + uux + -px = vuxx. (58)

P

Пусть стенка движется со скоростью U в положительном направлении оси Ox. Введём систему координат

t

x' = x — J Udt, (59)

о

движущуюся вместе с границей, и перепишем в ней уравнение (58):

ut + (u - U)ux' + 1 Px' = vux'x'. (60)

P

Отнесём это уравнение к движущейся стенке. В силу условия прилипания (57) на стенке u = U получим условие для давления в самом общем виде

Px' = P (-Ut + vux'x') . (61)

Производная ux'x' рассчитывается от скорости жидкости на стенке по нормали к ней. При U = const и течении при больших числах Рейнольдса из (61) следует приближение пограничного слоя px' = 0, которое совпадает с условием для давления и в случае неподвижных сеток. Поэтому в настоящей работе для давления на подвижной твёрдой границе используется также приближение пограничного слоя.

4. Стандартная k — e-модель турбулентности на подвижных сетках

Выпишем модификацию метода решения уравнений k — е-модели [1] на подвижные сетки.

Каждое из уравнений k — е-модели может быть записано в общем виде

I+£ —(-+*) Ц)=^ («2)

Интегрируя соотношение (62) по подвижному объёму и применяя формулу

/ dtdV = dt S ^dV — / ^ ■ dS, (63)

V (t) V (t) dV

получим уравнение в виде интегрального закона сохранения

d j ydV = — У у (w — xt) ■ dS + У (u + VydS + J HvdV. (64)

V(t) dV (t) dV (t) ^ V (t)

Неявная аппроксимация интегрального уравнения (64) приводит к системе нелинейных уравнений

3 (yV)-+! — 4 Шj + (yVГ"1

(У )jk-, )jk + (У )jk = RHS^1, (65)

2At J

где значения ф отнесены к центру ячейки и

ИШ^1 = - £ [(фт • Б)т+1/2 - (^ • Б)т-1/2^ +

п+1

+ Е

+ X] ' Б)т+1/2 - (фХ* • Б)т-1/2] +

V \( Вф ВФ дф

'Л--+ $у— + Бг —

а^) \ дх ду дг

„ П+1

а дФ . с дФ . с

V +--Б^ — + Бу— + Бг —

а^) \ дх ду дг

т+1/2

+ (Н V П1. (66)

П+1

т-1/2

Схема для вычисления невязкого потока приобретает вид

Гт+1/2 = 1 [(Мт + (^+1) • Бт+1/2-

- (фт + фт+1) • Цдт+1/2 -|и - ид |т+1/2 (фт+1 - фт) , (67)

где и = w • Б, ид = х* • Б. Скорость движения грани ячейки идт+1/2 находится по формулам (41), (45).

5. Результаты расчётов

5.1. Расчёт однородного потока на подвижной сетке

Целью численного эксперимента, описанного в настоящем разделе, является подтверждение точного выполнения дискретного УГК в предложенном численном алгоритме. Для этого рассматривается однородный поток несжимаемой жидкости с единичной скоростью через куб, параллельный четырём его граням. Поскольку компоненты скорости потока и давление в нём постоянны, то погрешность аппроксимации разностной задачи будет определяться только погрешностью выполнения УГК (24). При этом не важен размер ячеек сетки. Необходимо продемонстрировать, что на существенно неравномерной сетке (рис. 2) погрешность аппроксимации остается равной нулю. Сетка имела

Рис. 2. Деформация сетки в кубической области течения однородного потока, проводимая на каждом шаге по £

Результаты расчёта однородного потока на подвижной сетке

Предложенный метод

Время с Исходный метод [1] с несогласованным расчётом с согласованным расчётом

0.0 1.0000000 1.0000000 1.0000000

0.4 0.8785331 1.0025764 1.0000000

0.8 0.7646654 1.0056505 1.0000000

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

Решение находилось двумя методами: исходным [1], не учитывающим движение сетки, и предложенным методом (см. раздел 2.5) с несогласованным и согласованным расчётом объёмов Уп и УП±1/2. С помощью каждого метода было осуществлено пять шагов по времени с ДЬ = 0.2 с. При этом на каждом новом шаге проводилась дополнительная деформация сетки. В таблице представлены рассчитанные скорости потока для каждого из применённых методов в разные моменты времени.

5.2. Движение кругового цилиндра в покоящемся однородном поле несжимаемой вязкой жидкости

В данном разделе исследуются свойства и возможности построенного метода расчёта течений на подвижных сетках на задаче ламинарного обтекания кругового цилиндра, решение которой содержит большой набор достаточно сложных гидродинамических явлений. Адекватное предсказание последних является важным требованием, предъявляемым к численным методам. Особенность проводимых исследований — нестационарная постановка задачи, в которой в отличие от большинства используемых не жидкость обтекает неподвижный цилиндр, а цилиндр движется по покоящейся жидкости. Отметим, что важность этого теста состоит в демонстрации эффективности и надёжности созданного программного инструментария.

5.2.1. Постановка задачи

Движение цилиндра диаметра d по области с покоящейся жидкостью задаётся посредством перемещения со временем части границы расчётной области, совпадающей с контуром цилиндра. Для построения сетки используется криволинейная система координат, нормально связанная с поверхностью цилиндра. Таким образом, вместе с цилиндром движется и сетка (рис. 3). В качестве неподвижной внешней границы расчётной области принимается окружность радиуса Я с центром в начале декартовой системы координат х, у. Цилиндр движется вдоль оси Ох справа налево. В момент Ь = 0 центр цилиндра находится в точке (0.8Я, 0), в конечный момент времени Ь = Т — в точке (—0.8Я, 0). Движение начинается мгновенно с постоянной скоростью и.

Пусть 0 < £ < 1 есть продольная, а 0 < п < 1 — поперечная координаты криволинейной системы координат, нормально связанной с цилиндром. Связь между декартовой и криволинейной системами координат задаётся формулами

Г х (£, п, Ь) = —г (£, п, ¿) с°8 2п£ + х0 (Ь) I У (£,П,Ь) = г (£,n,t)sin2п£,

Рис. 3. Схема движения кругового цилиндра в области покоящейся несжимаемой вязкой жидкости: а — начальное, б — конечное положения цилиндра

где xo (t) = 0.8R — Ut — координата двигающегося по оси Ox центра цилиндра;

Ч d fw. ч d\ (1 + a)n — 1

L (f, t) = yjR2 — (xo (t) sin2nf)2 + xo (t) cos 2nf.

Кроме рассмотренных выше краевых условий на подвижной поверхности цилиндра, задаются все параметры невозмущённого потока на внешней границе расчётной области. На совпадающих границах Г+ и Г- (см. рис. 3) ставятся условия периодичности.

5.2.2. Физические и схемные параметры задачи

Для расчёта приняты следующие значения физических параметров задачи: d = 1 м, и = 0.125 м/с, коэффициент кинематической вязкости V = 0.003125 м2/с, давление в покоящейся жидкости Р = 0. Этому режиму обтекания цилиндра соответствовало число Рейнольдса

Ud л Яе = — = 40.

V

К схемным параметрам относятся радиус внешней границы Я = 200 м, время движения цилиндра Т = 2560 с, параметр сгущения сетки по нормальному к цилиндру направлению а = 50, шаг по физическому времени Д£ = 0.1 с, коэффициент искусственной сжимаемости в = 4, шаг по псевдовремени Дт = 1 с.

Сетка имеет 500 узлов в окружном и 500 узлов в нормальном к цилиндру направлениях. Заданное значение параметра сгущения а обеспечивало отношение нормального размера ячейки у внешней границы к размеру ячейки у поверхности цилиндра, равное 50.

а

5.2.3. Результаты расчёта

Реализуемому в расчёте числу Рейнольдса Яе = 40 соответствует стационарный режим ламинарного обтекания цилиндра. В то же время расчёт проводится в рамках неста-

ционарной постановки задачи, в которой на каждом слое по времени Ь решение устанавливается по псевдовремени т. Критерием сходимости итераций з по псевдовремени было условие

И

где норма определяется как

*)П+^ - 4 ^)п + ^)п-1

2АЬ

(ИИБ *)

п+1,в

< 10

-в

(69)

Е

тах ^|.

Здесь т = 1,..., 4 — номер уравнения, 1,3,к — номера ячеек по соответствующим координатным направлениям.

При выполнении условия (69) осуществляется переход на следующий слой по физическому времени и продвижение цилиндра по оси Ох влево. В процессе движения цилиндра в его окрестности формируется и устанавливается поле течения, характерное для обтекания неподвижного цилиндра потоком жидкости с соответствующими параметрами. По достижении цилиндром центра расчётной области решение становится стационарным и картина обтекания цилиндра приобретает вид, изображённый на рис. 4, где хорошо видны две сформировавшиеся симметричные рециркуляционные зоны на подветренной стороне цилиндра.

Полярный угол 9 точки отрыва потока с цилиндра и длина рециркуляционной зоны Ь в зависимости от положения цилиндра показаны на рис. 5. Коэффициент сопротивления

2 /V 2 ди

С Б = — I Р +

ри2

Ие дп

пх dБ

и коэффициент давления в передней критической точке

2(р - Р)

Ср =

ри2

показаны в процессе движения цилиндра на рис. 6. Здесь же серыми полосами отмечены интервалы изменения указанных параметров в сводных данных [14-16], а штрихом — значения этих параметров, рассчитанные в постановке "неподвижный цилиндр в потоке жидкости" в стационарном приближении на неподвижной сетке [1]. Видно хорошее совпадение параметров потока, полученных в нестационарной постановке на движущейся сетке, с экспериментальными данными и с результатами расчёта стационарного обтекания цилиндра на неподвижной сетке.

Рис. 4. Картина установившегося обтекания цилиндра (линии тока): 9 — угол отрыва потока, Ь — длина рециркуляционной зоны

о

Х0Ш

Рис. 5. Зависимости в (а) и Ь (б) от положения центра цилиндра: серые полосы — интервалы изменения значений параметров в сводных данных [14-16], штрих — расчёт в стационарной постановке [1]

Рис. 6. Зависимости С в (а) и Ср (б) от положения центра цилиндра: серые полосы — интервалы изменения значений параметров в сводных данных [14-16], штрих — расчёт в стационарной постановке [1]

Необходимо отметить, что рассмотренный подход при решении данной задачи является более затратным по времени по сравнению с традиционным. Кроме того, как следует из рис. 5, 6, движение цилиндра сопровождается осцилляциями параметров рассчитываемого потока, природа которых связана с изменением на каждом шаге по времени аппроксимирующей поверхность цилиндра ломаной кривой — многогранника. Верши-

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

5.3. Моделирование течения в гидротурбине

при закрытии лопаток направляющего аппарата

Рассматривается задача, для решения которой и создавался предложенный здесь метод. Эта задача в принципе не может быть решена с использованием подходов, базирующихся на фиксированных сетках. Отметим преимущества представленного метода при решении данной задачи. Во-первых, в методе сохранены все примечательные свойства исходного численного алгоритма моделирования течений на фиксированных сетках [1]. Для высокой экономичности численного решения трёхмерных нестационарных уравнений Рейнольдса несжимаемой жидкости в реальных аэрогидродинамических установках, имеющих сложные многосвязные проточные части, в алгоритме использован максимально устойчивый и в то же время достаточно просто реализуемый неявный метод конечных объёмов. Экономичность алгоритма повышается тем, что сегменты, в которых находится решение, строятся топологически эквивалентными параллелепипедам с регулярными сетками в них. Благодаря этому, а также введению фиктивных слоев сетки алгоритмы остаются однородными во всем сегменте, в том числе и у его границ. Обращение неявных операторов посредством попеременно-треугольного метода становится в данном случае крайне простым. Сложная область проточного тракта с помощью созданного аппарата геометрического моделирования покрывается указанными сегментами. При этом разработан оригинальный экономичный метод задания и хранения в памяти компьютера типов граничных условий для каждого сегмента. Метод позволяет задавать произвольные граничные условия на различных участках границы и организовывать оперативный обмен данными между сегментами при минимальных затратах памяти и времени выполнения операций. Проводимая естественным образом сегментация проточного тракта аэрогидродинамической установки позволяет, с одной стороны, эффективно находить решения уравнений в каждом сегменте, с другой — строить в каждом сегменте достаточно качественные сетки, независимые от сеток соседних сегментов.

Помимо сохранения свойств исходного алгоритма, метод дополнен очень важным для работы на подвижных сетках свойством точного выполнения УГК.

Решается задача моделирования нестационарных переходных процессов в гидравлической турбине. Переходной процесс — смена одного установившегося течения в гидротурбине другим, вызванная изменением угла установки лопаток направляющего аппарата. При этом меняется расход воды, пропускаемый турбиной, что приводит к появлению гидравлического удара в длинном подводящем водоводе и росту или уменьшению скорости вращения рабочего колеса. Необходимость учёта этих явлений делает задачу моделирования переходных процессов существенно более сложной, нежели моделирование отдельного установившегося режима работы.

Моделирование упругого гидроудара в подводящем водоводе проводится путём решения одномерных уравнений акустики для области водовода. Расчётная область для гидротурбины, в которой решаются нестационарные уравнения движения несжимаемой жидкости (11), состоит из направляющего аппарата, рабочего колеса и отсасывающей трубы (рис. 7). Способ постановки краевых условий на входе и выходе, при которых расход жидкости находится из решения задачи, а также способ обмена параметрами потока между водоводом и самой гидротурбиной подробно описаны в [17].

Ниже представлены результаты моделирования одного из переходных процессов — уменьшения мощности турбины. В начальный момент времени £ = 0 задаётся стационарное поле течения, полученное при расчёте режима максимального КПД. Далее, за 10 с от момента времени £ = 2.5 с, угол открытия лопаток НА уменьшается линейно от 26.6° (соответствует режиму максимального КПД) до 18.9° (режим неполной загрузки). Шаг по времени Д£ = 0.01 с соответствует повороту РК на 12°. Подвижная сетка в канале НА на каждом шаге по времени строилась алгебраически (рис. 8).

Известно, что в режиме максимального КПД поток за РК имеет слабую закрутку, вследствие чего течение в конусе отсасывающей трубы приобретает практически стаци-

Рис. 7. Расчётная область и сетки в направляющем аппарате, рабочем колесе и отсасывающей трубе гидротурбины

Рис. 8. Сетки в межлопаточном канале направляющего аппарата перед началом его закрытия, £ = 2.5 с (а) и после закрытия, £ = 12.5 с (б)

онарный осесимметричный характер. В режиме неполной загрузки, напротив, окружная скорость потока за РК существенна, что приводит к формированию прецессирую-щего вихревого жгута в конусе ОТ.

На рис. 9 представлено полученное в расчёте изменение момента сил, действующих на РК, и расхода во времени, на рис. 10 показана эволюция давления перед НА и на стенке конуса ОТ. Снижение расхода при закрытии лопаток приводит к положительному гидроудару, который проявляется в повышении давления перед НА (рис. 10, а). Зародышевый вихревой жгут, наблюдающийся в начальный момент времени, по мере закрытия НА увеличивает свой размер и радиус прецессии (рис. 11), что вызвано ростом остаточной закрутки за рабочим колесом. После окончания движения лопаток НА, т.е. при £ > 12.5 с, течение имеет периодически нестационарный характер. Вращение вихревого жгута вызывает пульсации расхода и момента сил на валу гид-

Рис. 9. Поведение момента сил М и расхода Q в процессе уменьшения мощности а б

Рис. 10. Эволюция давления перед НА (а) и на стенке в конусе ОТ (б)

= 6.02 с

¿2 = 8.06 с

¿з = 10.02 с

и = 12.80 с

Рис. 11. Рост интенсивности вихревого жгута (показаны изоповерхности давления в моменты времени ¿1 , ¿2, ¿3,¿4, отмеченные на рис. 10, б)

ротурбины (рис. 10, б). Амплитуда пульсаций давления в конусе ОТ достигает 4% от действующего напора. Период прецессии жгута Т = 1.41 с, что согласуется со значением, измеренным в эксперименте Техр = 1.3 с.

Заключение

Метод расчёта течений несжимаемой жидкости в неподвижных областях [1] обобщён на задачи с подвижными границами. Предложен подход, обеспечивающий точное выполнение условия геометрической консервативности на дискретном уровне. Условиями, накладываемыми на движущуюся сетку, являются неизменность числа ячеек сетки и ограниченность длины передвижения ячейки за один шаг по времени окрестностями этой ячейки (новое её положение должно иметь общую часть со старым положением). Построенный метод применён для решения тестовых и практических задач. Анализ результатов решений показал, что метод корректно воспроизводит все рассмотренные характеристики и может быть использован для моделирования течений в областях с подвижными границами в различных задачах вычислительной гидродинамики.

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

[1] Чёрный С.Г., Чирков Д.В., Лапин В.Н. и др. Численное моделирование течений в турбомашинах. Новосибирск: Наука, 2006. 202 с.

[2] Trülio J.G., Trigger K.R. Numerical Solution of the One-Dimensional Hydrodynamic Equations in an Arbitrary Time-Dependent Coordinate System. Techn. Rep. UCLR-6522. Univ. of California Lawrence Radiation Laboratory, 1961.

[3] Thomas P.D., Lombart C.K. Geometric conservation law and its application to flow computations on moving grids // AIAA J. 1979. Vol. 17, No. 10. P. 1030-1037.

[4] Demirdzic I., Peric M. Space conservation law in finite volume calculations of fluid flow // Intern. J. Numer. for Methods in Fluids. 1988. Vol. 8. P. 1037-1050.

[5] Shyy W., Udaykümar H.S., Rao M.M., Smith R.W. Computational Fluid Dynamics With Moving Boundaries. Washington, DC: Taylor and Francis, 1996.

[6] Lesoinne M., Farhat C. Geometric conservation laws for flow problems with moving boundaries and deformable meshes, and their impact on aerolastic computations // Comput. Methods in Appl. Mech. and Eng. 1996. Vol. 134. P. 71-90.

[7] Koobüs B., Farhat C. Second-order time-accurate and geometrically conservative implicit schemes for flow computations on unstructured dynamic meshes // Ibid. 1999. Vol. 170. P. 103-129.

[8] Forster C., Wall W.A., Ramm E. On the geometric conservation law in transient flow calculations on deforming domains // Intern. J. for Numer. Methods in Fluids. 2006. Vol. 50. P. 1369-1379.

[9] Engel M., Griebel M. Flow simulation on moving boundary-fitted grids and application to fluid-structure interaction problems // Ibid. 2006. Vol. 50. P. 437-468.

[10] Зайцев Н.О., Щур Н.А. Применение деформируемых сеток для численного моделирования течений в областях с подвижными границами // Научно-техн. ведомости Санкт-Петербургского гос. ун-та. 2006. Т. 47. С. 15-22.

[11] Волков К.Н. Дискретизация уравнений Навье — Стокса на подвижных неструктурированных сетках // Вычисл. методы и программирование. 2008. Т. 9. С. 256-273.

[12] Фихтенгольц Г.М. Курс дифференциального и интегрального исчисления. Т. 2. М.: Физматлит, 2001. 810 с.

[13] Anderson W.K., Thomas J.L., van Leer B. Comparison of finite volume flux vector splittings for the euler equations // AIAA J. 1986. Vol. 24, No. 9. P. 1453-1460.

[14] Takami H., Keller H.B. Steady two-dimensional viscous flow of an incompressible fluid past a circular cylinder // Phys. Fluids Suppl. 1969. Vol. 12. P. 11-51.

[15] Coütanceaü M., Bernard R. Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 1. Steady flow // J. of Fluid Mech. 1977. Vol. 79. P. 231-254.

[16] Tritton D.J. Experiments on the flow past a circular cylinder at low Reynolds numbers // Ibid. 1959. Vol. 6. P. 547-567.

[17] Cherny S.G., Chirkov D.V., Bannikov D.V. et al. 3D Numerical simulation of transient processes in hydraulic turbines // 25th IAHR Symp. on Hydraulic Machinery and Systems. Timisoara, Romania, 2010. P. 1-8.

Поступила в 'редакцию 25 ноября 2011 г., с доработки — 31 июля 2012 г.