ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПРОСТРАНСТВЕННЫХ ТУРБУЛЕНТНЫХ ТЕЧЕНИЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ НА ОСНОВЕ K - є МОДЕЛЕЙ Текст научной статьи по специальности «Физика»

ПАМЯТИ АКАДЕМИКА АНАТОЛИЯ ФЕДОРОВИЧА СИДОРОВА

(1933-1999)

31 марта 1999 года ушел из жизни крупнейший ученый-математик, член Президиума Уральского отделения РАН, директор Института математики и механики УрО РАН, академик РАН Анатолий Федорович Сидоров.

А. Ф. Сидоров в 1955 г. окончил Ленинградский государственный университет. В Институте математики и механики УрО РАН работал с апреля 1963 г. в должности заведующего отделом прикладных задач. В 1963 году А. Ф. Сидоров защитил кандидатскую, а в 1969 — докторскую диссертации. В 1971 г. ему было присвоено ученое звание профессора. В 1987 г. А. Ф. Сидорова избирают членом-корреспондентом АН СССР, а в 1991 г. действительным членом РАН.

А. Ф. Сидоров — крупный ученый в области механики и вычислительной математики. Его творчество в области математики было многогранно. Одно из основных направлений, в котором А. Ф. Сидоров достиг выдающихся результатов, — это аналитические методы исследования в газовой динамике и гидродинамике. Здесь Анатолий Федорович получил общие результаты в теории бегущих волн, впервые вывел уравнение тройных волн, построил серии точных решений, доказал теоремы о примыкании бегущих волн различных рангов.

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

Результаты научных разработок А. Ф. Сидорова широко используются российскими организациями в практике конструирования и создания объектов новой техники.

А. Ф. Сидоров был одним из инициаторов развития в России нового перспективного направления суперкомпьютеров МВС-100 — МВС-1000, включаю-

щего создание аппаратных и программных средств и разработку математических алгоритмов параллельного действия для решения задач, требующих больших вычислительных мощностей. В ИММ УрО РАН на базе МВС-100 создан современный информационно-вычислительный центр. Под руководством А. Ф. Сидорова в Уральском регионе развернута деятельность по телекоммуникационному обеспечению УрО РАН и выходу в ИНТЕРНЕТ.

А. Ф. Сидоровым опубликовано более 140 научных работ, включая одну монографию. Возглавляя отдел прикладных задач, А. Ф. Сидоров являлся одновременно руководителем нескольких крупных тем, выполняемых по постановлениям директивных органов и получивших высокую оценку.

Большую научную работу Анатолий Федорович успешно сочетал с педагогической деятельностью в качестве профессора Уральского государственного университета и заведующего кафедрой параллельных компьютерных технологий УрГУ при ИММ, проводил большую работу по подготовке высококвалифицированных кадров. В числе его учеников два доктора и 17 кандидатов наук.

Анатолий Федорович Сидоров являлся Председателем Объединенного Ученого Совета по математике, механике и информатике УрО РАН, членом Президиума УрО РАН, членом специализированных советов по защитам диссертаций при ИММ и УрГУ, членом редколлегии журналов "ЖВМ и МФ", "Вычислительные технологии", "Квант", "Theoretical and Applied Mechanics", "Numerical Analysis and Mathematical Modeling". В 1983 г. он был избран членом Национального комитета по теоретической и прикладной механике. В последние годы А. Ф. Сидоров являлся членом Совета Российского фонда фундаментальных исследований.

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

Редакционная коллегия

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

Том 4, № 2, 1999

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПРОСТРАНСТВЕННЫХ ТУРБУЛЕНТНЫХ ТЕЧЕНИЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ НА ОСНОВЕ k - б МОДЕЛЕЙ *

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

An efficient numerical method is considered for calculating three-dimensional turbulent flows of an incompressible fluid based on the artificial compressibility method for the solution of the Reynolds equation. For the closure of the latter various modifications of a two-parameter k — e turbulence model are employed. The equations are approximated by the implicit finite volumes method and the third order TVD-schemes. The constructed schemes are implemented by the alternately triangular method. The results of the numerical calculations are presented.

Введение

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

Целью настоящей работы является распространение предложенного в [1] метода на уравнения Рейнольдса с замыкающей двухпараметрической k — б моделью турбулентности. Данная модель относительно проста, что играет немаловажную роль при моделировании сложных пространственных турбулентных течений. В то же время незначительная модификация модели, например, учет анизотропности тензора напряжений подобно тому, как это делается в релаксационной k — б модели Дурбина [2], или внесение поправочного на вращение источникового члена [3] даст альтернативный сложным моделям подход, не уступающий им по точности.

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, гранты №98-01-00742 и 96-01-00136, и РФФИ-ИНТАС, грант №95-1149.

© С.Г. Черный, П.А. Шашкин, Ю.А. Грязин, 1999.

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

Используемая к — е модель турбулентности справедлива для полностью развитых турбулентных течений. Вблизи же твердой поверхности турбулентные флуктуации подавляются стенкой и течению присущи как ламинарные, так и турбулентные свойства. Известная ее модификация для низких чисел Рейнольдса, решающая эту проблему посредством демпфирующих функций [4], неприемлема, так как для ее использования требуется большое число узлов в узком слое возле твердой стенки, что в трехмерном случае приводит к огромным затратам ресурсов компьютера. Поэтому за основу берется И^О к — е модель турбулентности [5], дополненная эмпирическим законом о поведении потока вблизи твердых границ (метод пристеночных функций).

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

1. Математическая модель 1.1. Уравнения Рейнольдса

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

Ъ = 0 (1)

dui d д * д

Ht + dXj (uiuj) + dx;p = dXj ) + fi' (2)

где

* 2. dui duj

P = P k ТИ =ô--+ , Ve = v + Vt,

3 dxj dxi

k — кинетическая энергия турбулентности, vt — турбулентная вязкость, p — осредненное давление, ui — осредненные компоненты вектора скорости в декартовой системе координат. Для расчета течений во вращающихся областях (например, в рабочем колесе гидротурбины), а также в поле тяжести в уравнения добавлен вектор массовых сил:

f = (fi, f2, f3) = —w x w x r — 2w x u + g,

включающий в себя центробежные, кориолисовы и гравитационные члены. Здесь w — вектор угловой скорости вращения.

1.2 Замыкающая модель турбулентности

Для вычисления турбулентной вязкости рассматриваются два варианта двухпараметри-ческой k — е модели турбулентности.

1.2.1. RNG k — е модель турбулентности

RNG модель турбулентности (renormalization group based k — е model) [5] получила в последнее время широкое распространение благодаря хорошему совпадению получаемых численных результатов с имеющимися экспериментальными данными, а также высокой скорости сходимости базового алгоритма. В этой модели параметры турбулентности вычисляются из уравнений

dk д

дк

де д

де

k2

-Trr^^—ikkUj - vk—— ) = Hk, —+ —(еи, - ) = He, Vt = CM —,

дЬ дх,

"дх

дЬ дх

дх,-

(3)

где е — скорость диссипации турбулентности,

Hk = G — е, He = С*!7 G — Се2~Г

kk

,дщ , ди, ди G = vXT- + 7Г2)

дх, дхг дхj

С* = С П(1 - П/По) п

Се1 = Се1---, , 0..3 , П

Vt

Vk = V +--, Vf

Jk

1 + вп3 Vt

G

С,е

0.5

V +

of

Эмпирические константы в приведенных уравнениях равны:

С = 0.0845, Се 1 = 1.42, Са = 1-68,

ак = 0.72, ае = 0.72, по = 4.38, в = 0.015.

1.2.2. Метод пристеночных функций для турбулентного пограничного слоя около

стенок в RNG модели

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

Ближняя к стенке зона, в которой вязкие напряжения доминируют над рейнольдсо-выми, называется вязким, или ламинарным, подслоем. Здесь имеет место линейная зависимость тангенциальной составляющей скорости от расстояния до стенки по нормали:

2

и+ = y+.

(4)

В (4) введены обозначения для безразмерной тангенциальной составляющей скорости

и+ = ^ (5) ит

и безразмерного расстояния до стенки

у+ = у— • (6)

V

Модуль тангенциальной составляющей скорости qт определяется по формуле

Ят =| Чт | = | и - (и ■ п-)п- (7)

где п— — единичный внешний нормальный к стенке вектор, ит — динамическая скорость, определяемая как

ит = л / Т„,,

где — напряжение трения на обтекаемой поверхности. Зависимость (4) эквивалентна по своей сути определению по линейной аппроксимации профиля скорости в слое между стенкой и пристеночным узлом Р (соответствует ламинарному режиму течения):

т*ш = V—. (8)

Ур

Вязкий подслой с линейным распределением скорости (4) в практике расчетов обычно определяется в интервале

0 < у+ < 11.6.

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

и+ = -1- 1п(Еу+) (9)

каг

и рейнольдсовы напряжения намного превышают вязкие. Здесь постоянная Кармана каг принимает значения

каг = 0.4 ^ 0.42.

Постоянная, определяющая степень шероховатости стенки, может изменяться в пределах

Е = 8.8 ^ 9.79.

Логарифмический подслой турбулентного пограничного слоя обычно определяется в интервале

11.6 < у+ < 400.

Область с у+ > 400 считается режимом полностью развитого турбулентного течения.

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

и2

кр = , (10)

и3

бр каг ур' (11)

Для определения в (10) и (11) динамической скорости ит перепишем (9) следующим образом:

а 1

£ = каг1п(Е Ке >■ (12)

Обозначив в (12) искомую величину динамической скорости ит через х, запишем уравнение для ее определения:

Р(х) = -1- 1п(£И,е у„х) - — = 0. (13)

каг х

Уравнение (13) решается с использованием итерационного метода Ньютона:

Р (хт)

__, ^т)

хт+1 хт

Р '(Хт)

где

г

^'(Х) = ка^ + ^ каг х х2

т — номер итерации.

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

ТГ) = Г*-. (14)

ду/ р каг ур

Выполнение последнего равенства гарантирует правильный расчет напряжения трения на стенке.

Другой путь реализации метода пристеночных функций состоит в объединении законов (8) и (9) в одно универсальное выражение

т- = V- —, (15)

Ур

где

каг У+

= V тах

1п(Еу+)_

Безразмерное расстояние у+ определяется при этом по формуле

С 1/4к1/2У

у+ = V Ур ■ (16)

Обобщенный коэффициент кинематической вязкости vш используется для подсчета напряжения трения на стенке.

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

1.2.3. к — е модель турбулентности с демпфирующими функциями для низких чисел

Рейнольдса

Данный вариант к — е модели представляется в виде

^ к2 дк д / дк\ ТТ де д / де\ тт

щ = /Лт, ж + дх~ - ^д; = Як' д* + дх; - *д; = Я- (17)

где

к е е2 е

Як = С - е - 2^—, Я = С£1 -С - С/т - /3.

у2 к к у2

Здесь у — кратчайшее расстояние до ближайшей твердой стенки. Демпфирующие функции определяются следующим образом [4]:

к2

/ = ехр [-2.5/(1 + Д*/50)], /2 = 1 - 0.3ехр (-Щ2), /з = ехр (-0.5у+), Щ = -.

Константы принимаются равными:

С = 0.09, СС1 = 1.44 ^ 1.59, Се2 = 1.9 ^ 2, ак = 1, ае = 1.3 ^ 1.47.

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

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

2.1. Метод решения уравнений Рейнольдса

Для дискретизации исходных уравнений (1), (2) используется неявный метод конечных объемов [1]. Вычисление потоков на гранях расчетной ячейки осуществляется таким образом, что результирующая схема имеет противопоточную аппроксимацию третьего порядка для конвективных членов и центрально-разностную второго порядка для вязких членов. Линеаризация полученных соотношений проводится с помощью метода Ньютона, причем в неявном операторе оставляются только члены, отвечающие за первый порядок аппроксимации невязких потоков и соответствующие повторным производным вязкой части потока. Решение полученной линеаризованной системы разностных уравнений осуществляется с помощью оригинального варианта попеременно-треугольного метода (ЬИ-факторизации). Отличительной чертой данного алгоритма является его абсолютная устойчивость и экономичность. Подробно алгоритм описан в работе [1].

2.2Метод решения уравнений для турбулентной кинетической энергии и скорости ее диссипации

2.2.1. Обобщенная запись замыкающих уравнений Каждое из уравнений к — е модели может быть записано в следующем общем виде:

дф + д (Ф * дф где ф, V* и Я приведены в таблице.

Яг + - V*£т) = Я, (18)

Таблица

Уравнение Ф v * Н

Для турбулент- к v + VI / ак С - е ЯУС модель

ной кинетичес- С - е - 2vk/y2 Модель для

кой энергии низких чисел Рейнольдса

Для скорости е V + VI/а* С*1 • е/k • с - С*2 • е2/k ЯУС модель

диссипации турбулентной

кинетической Се1 • е/k • С - С*2 • е2/k • ¡2 - 2v • е/у2 • / Ыодель для

энергии низких чисел Рейнольдса

2.2.2. Дискретизация уравнений

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

I Ф^У = - £ фи^Б + ^ V^ф^Б + У Я^У, (19)

V дУ дУ V

где дУ — замкнутая поверхность произвольного фиксированного объема V; ^Б = п • — элемент поверхности дУ, умноженный на единичную внешнюю нормаль п к ней.

Рис. 1. Текущая расчетная ячейка.

Расчетная область разбивается на элементарные ячейки в виде криволинейных шестигранников (рис. 1). Неявная аппроксимация интегрального уравнения (19) на текущей

расчетной ячейке приводит к следующей системе нелинейных уравнений

дга+1 фп

"л*

Pijk , . PljkVj = (20)

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

1 = - £ [(фи ■ S)m+1/2 - (фи ■ S)m-1/2]n+1 +

m=i,j,k

+ ^ ([v*(S* + Sy ^ + )]m+1/2 - [v*(S* + Sy ^ + )]m—1/2 Г+1 + Hij+ Vj;

Sx, Sy, Sz — компоненты вектора Sm+1/2 в декартовой системе координат.

2.2.3. Определение разностных потоков. Невязкие потоки

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

(фи ■ S)m+1/2 = 1[((фи)т + (фи)т+1) ■ Sm+1/2- | U |т+1/2 Дт+1/2ф] - Wm+1/2, (21)

U = и ■ S, | U |= U+ - UДт+1/2ф = фт+1 - фт,

U+ = 0.5(U + | U |), U- = 0.5(U-| U |).

В (21) член в квадратных скобках аппроксимирует невязкий поток на грани m + 1/2 c первым порядком; Wm+1/2 — член, повышающий порядок аппроксимации невязкой части схемы (20) до второго или третьего:

1 — Ф 1 + Ф

Wm+1/2 = —[N-+3/2 - M+-1/2] + — M-+1/2 - Nm++1/2],

Nm+3/2 = minmod(U-+3/2^+3/2^eU-+1/2^+1/2^, m—1/2 -

minmod(U+ —1/2 Д m— 1/2^^+1/2 Дт+1/2ф) ,

Mm+1/2 = minmod(Um+1/2^+1/2^eUm+3/2 Дm+3/2ф),

m+1/2

— 1/2Дт— 1/2ф) ,

minmod(x,y) = sign(x) max[0,min(| x |,y ■ sign(x))],

1 < в < (3 - Ф)/(1 - Ф), m = i,j,fc. (22)

Без минмодульных ограничителей схема (20) - (22) имеет для невязких членов направленные против потоков разности второго (Ф = -1) или третьего (Ф = 1/3) порядка аппроксимации.

Вязкие потоки

Для аппроксимации на грани ячейки т +1/2 (т = г, к) вязкого потока

5ф 5ф

5ф,

К дХ + °у дУ + 5 ^)]т+1/2

введем условную криволинейную систему координат

С = V = С = С

(23)

(24)

связанную с разбиением расчетной области на элементарные ячейки. Координатные линии С соответствуют изменению индекса г при фиксированных ] и к; линии п соответствуют изменению индекса ] при фиксированных г и к; линии £ соответствуют изменению индекса к при фиксированных г и ^ (см. рис. 1). Далее представим все производные, входящие в (23), через производные по введенным криволинейным координатам:

5 5С 5 + + 5( 5 5С '

5 5С 5 + 5 + 5( 5 5у 5у 5С 5('

А = 5СА + + (25)

Вместе с тем можно найти связь между коэффициентами преобразования производных и геометрическими параметрами сетки:

/ 5С 5С 5С \ =

\5ж' 5у' 5г/ г+1/2 /2

/ 5т Дп

¿+1/2 ^,+1/2

/ 5С 5С 5С \ ДС

(5х ) 5у ) )г+1/2) (°Х) 5у) )¿+1/2;

(°Х) °у, )й+1/2;

(26)

ЧЗя'Зу'Зг/к+1/2 Н+1/2

где ^+1/2 = 0.5(Ут + ^+1).

Переходя в (23) от производных по декартовым координатам к производным по криволинейным координатам по формулам (25), получим

где

[V '(5, ^ 1 4 + + ^ 5Ф )" , = J т+1/2 И- |

= ^ °х + 5С 5 + 5С

^2 = ^с + О °Х + + 5у у

= + ^ °Х + 5С 5 + „ Оу +

т+1/2

(27)

5, 5.

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

Грань 1+1/2

Аппроксимация производных:

дф , ,

~ фг+1;,к - фг,;,к,

дф

— ~ 0.25(фг ;+1,к - фг; —1,к + фг+1;+1,к - фг+1; —1,к),

дф дс

Вычисление коэффициентов преобразования:

де

дд

дп дд

дс

дд

где д есть любая из декартовых координат х, у и г

0.25(фг,;,к+1 - фг;,к—1 + фг+1,;,к+1 - фг+1,;,к-1).

— ~ )г+1/2,;,к/^г+1/2,;,к,

0.25[(5,)г,;+1/2,к + )г+1,;+1/2,к + )г,;-1/2,к + )г+1; —1/2,к] / ^¿+1/2,;,к,

0.25 [(5,) г,;,к-1/2 + )г+1,;,к-1/^ / ^¿+1/2,;,к,

Грань 2+1/2

Аппроксимация производных:

дф

- 0.25(фг+1,;,к - Фг-1,;,к + фг+1;+1,к - фг—1;+1,к), дФ^ф л

дП ~ фг,;+1,к - фг,;,к,

— « 0.25(фг,;,к+1 - фг,;,к—1 + фг;+1,к+1 - фг,;+1,к-1).

5ф дС

Вычисление коэффициентов преобразования:

— ~ 0.25[(5',)г+1/2,;,к + )г+1/2,;+1,к + )г— 1/2;,к + )г— 1/2;+1,к] / ^¿;+1/2,к,

дп дд

дС

дд

0.25[(5,)г,;,к+1/2 + )г,;+1,к+1/2 + )г,;,к—1/2 + )г,;+1,к—1/^ / ^;+1/2,к.

Грань k+1/2

Аппроксимация производных:

дф

дф

— ^ 0.25^i,j+1,fc — фг— i+ фг,7-+1,к+1 — фг,7 —i,fe+i),

ф

Вычисление коэффициентов преобразования: 0.25 [(Sq)

г+1/2 + (Sq )i

+ 1/2 j,k+1

+ (Sq )i

дп дд

дС

дд

0.25 [(Sq ) + (Sq ) + (Sq ) + (Sq ) ij —1/2,k+iJ / Vi,j,k+1/2>

— ~ (Sq)i,j,fc+i/2/Vi,j,fc+i/2.

Другой способ аппроксимации на грани ячейки m +1/2 производных дф/дх, дф/ду и дф/ дг, входящих в выражение для вязкого потока (23), заключается во введении контрольного объема Vm+i/2, содержащего внутри себя данную грань. Например, на рис. 2 сечение плоскостью k = const основного объема Vjk обведено жирной линией, а контрольного объема Vi+i/2!j,fe возле грани i + 1/2 — заштриховано. Тогда, используя следствие формулы Остроградского — Гаусса (см. выражения (30) в разделе 2.2.4), получим для грани i + 1/2:

(дИ г+i/2 = (фSq )i+i'j'k — (фSq )ij>fc+

+ (фSq ) г+1/2 j+1/2,k — (ФSq ) г+1/2 j —1/2,k + ^Sq ) i+1/2 j,k+1/2 (ФSq ) i+1/2,j,k—1/2.

• • •

i-l,j+l ij+l i+lj+l

• Щ

i-lj

• • •

ij-l i+lj-1

Рис. 2. Введение вспомогательного контрольного объема V^+i/2,j,fc•

2.2.4■ Аппроксимация источникового члена Н

Рассмотрим аппроксимацию источниковых членов Нд и Не в правых частях к — е модели (3) или (17).

Член генерации диссипации В выражения для Нд и Не входит член генерации диссипации:

2 / Ячп\ 2

С = 2

/дм\2 /дм\2 /ди>у \дж/ \ду/ V дг /

+

/дм дм \2 /дм дм> \2 /дм дм> у \ду дж/ \дг дж/ \дг ду/

(28)

В схеме (20) он аппроксимируется в центре ячейки г, ^, к на п-м слое по времени. Поэтому для вычисления производных от компонент скорости по декартовым координатам, входящих в С, требуются соотношения, которые будут отличаться от применяемых в разделе 2.2.3 для расчета производных от ф на гранях ячейки.

Для вычисления производных от компонент скорости, входящих в С, воспользуемся формулой Остроградского — Гаусса:

V ■ ф^У

Отсюда, последовательно полагая, что

J = £ ^¿Б.

V V дУ

(29)

получим соответственно

0, 0),

ф =4 (о,^, о), (о, 0,ф),

У дХЛ V = £ фйБх,

v ду

/ду^=/ ^,

v ду

дф дг

¿V = ф .

(30)

v ду

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

(ЦУ) = )^+1/2 - )г-1/2 + Ь+1/2 - ),-1/2 + Ь+1/2 - )к-1/2. (31)

Значения компонент скорости на гранях ячейки вычисляются по формулам

Фт+1/2 = 0.5(^т+1 + фт). (32)

С учетом выполнения всегда интегрального равенства

¿Б = 0 (33)

дУ

и, как следствие, его дискретного аналога

§¿+1/2 — 8г-1/2 + 8^+1/2 — 8^-1/2 + Бд+1/2 — БД-1/2 = 0, (34)

а также соотношений (32), выражения для производных от компонент вектора скорости (31) могут быть упрощены:

( д"V ) = 0.5(^+15^+1/2 - фг-1^91-1/2 + ^^+1^9^+1/2 - 'ф7-15ад-1/2 +

+фк+15дк+1/2 — фк-15дк-1/2). (35)

Явно-неявная аппроксимация источникового члена

Свободные члены в правых частях уравнений (3) или (17) к — е моделей имеют в общем случае вид

к

Н = С - е - 2^—,

V2'

е е2 е

Не = С*1 -С - Се2т/2 - 2^/^—^. (36)

к к у2

Для того, чтобы увеличить запас устойчивости численного алгоритма (20) решения к — е уравнений (3) или (17) (а следовательно ускорить сходимость численного решения этих уравнений к стационарному решению посредством выбора любого сколь угодно большого шага Д£), будем все входящие сюда слагаемые, имеющие отрицательные коэффициенты при искомых функциях, аппроксимировать неявно (т.е. на п + 1-м слое по времени). Перепишем выражения для Нд и Не в (36) следующим образом:

1к Н = С - тк - 2^-,

Н = С!-С - Се2/ - 2/з4, (37)

к т V2

где Т = к/е есть турбулентный масштаб времени. Обычно он ограничивается снизу значением

Т = тах

где

(38)

Су —

Поэтому коэффициенты - (1/Т + 2^/у2) при к в Нк и - (Се2/Т + 2^/у2/3) при е в Не являются всегда отрицательными. В силу этого только данные слагаемые в свободных членах аппроксимируются неявно, а остальные — явно:

1 \ п

1 V

НП+1 = С т + 2у^1 кп+1,

ЯГ1 = - ЦСТ2 + 2-V!/з) е

п+1

(39)

Следует отметить, что в свободных членах (39), как и в остальных членах разностных уравнений (20), составляющие вектора скорости берутся уже посчитанными на п + 1-м слое из уравнений неразрывности и количества движения.

2.2.5. Линеаризация

Для решения нелинейного уравнения (20) проводится его линеаризация по методу Ньютона с использованием записи линеаризованных уравнений в дельта-форме (т. е. уравнение записывается относительно Д^га+1 = фга+1 — фп):

— (

Д£ Удф )

(фп+1 — фп) =

(40)

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

(фи ■ Б)т+1/2 = 2 {[(Фи)т + (фи)т+1] ' Бт+1/2 — | и |т+1/2 Дт+1/2ф} ; (41)

а в вязком потоке оставляются производные только по той координате, поверхностью уровня которой является рассматриваемая грань:

" дф 06 дф

" дё + дП + дс.

' 06 06 дф"

" дё + дП + дс.

дф дф дф

" (£* дё + дП + дс.

¿+1/2

* дф дё

¿+1/2

¿+1/2

щ),

¿+1/2

* дф\ °сл

. (42)

к+1/2 V / к+1/2

Это означает, что, вообще говоря, в разностных потоках уравнений (20) неявно рассмотрены только члены, указанные в правых частях равенств (41) и (42), а остальные берутся явно.

Учитывая (41), (42), а также аппроксимацию свободных членов (39), найдем дЛЯ£/дф и перепишем (40) в виде

'V _ +

— + 0У + ит/2Дг+1/2 + и^_1/2Дг_1/2 — [(^^1) ¿+1/2 Д ¿+1/2 — )г_1/2Дг_1/2] +

+ и^_+1/2Д^+1/2 + и+-1/2Д^_1/2 — [К ^2) ¿+1/2 Д ¿+1/2 — О* ^ ¿-1/2 Д_1/2] + и_+1/2 Дк+1/2 + + и+_1/2 Дк_1/2 — [(V * ^3)^+1/2 Дк+1/2 — (V * из)*_1/2 Д*_1/2] } Д^ = , (43)

где

^ \ 1/Т + 2v/y2 для к-уравнения,

Се2 /Т + 2v/y2/3 для е-уравнения.

Далее идентифицируем коэффициенты (т = 1, 2, 3) так же, как и величины и±, только по нижнему индексу:

п

^¿±1/2 = (^1 )г±1 /2 — коэффициент ш>1, взятый на грани 5г±1/2; ^•±1/2 = (ш2).±1/2 — коэффициент ш2, взятый на грани 5^±1/2; ^^±1/2 = (^з)к±1/2 — коэффициент ш3, взятый на грани 5^±1/2.

Тогда уравнение (43) можно переписать в более компактном виде:

V - + - +

(А! + + Сг+1/2Дг+1/2 + Сг-1/2 Дг-1/2 + С.+1/2 Д.?+1/2 + С?-1/2Д^-1

+ Сд+1/2 Дк+1/2 + С+-1/2Дк-1/2)ДфП+1 = ДН5П, (44)

где

Ст+1/2 = и?±+1 /2 ± ^*ш)т+1/2, т = г, к.

Из определения шт в (27), а также с учетом (26) и введенной выше идентификации ш по нижнему дробному индексу получим выражение для шт+1/2 в (44):

8т+1/2 ■ 8т+1/2 ,

шт+1/2 = -^-. (45)

Угп+1 /2

2.2.6. Метод ЬЦ-факторизации решения линеаризованного уравнения Перепишем систему уравнений (44) следующим образом:

(^-С-1/2ТГ +С-+1/2Т+ -.-1/2. + С"+1/2Т?+ -С+-1/2Тк- + С-+1/2Т+) ДфП+1 = ^^ (46)

где V

П = + + С+1/2 - Сг+1/2 + С+-1/2 - С.+1/2 + С^-1/2 - Ск+1/2;

т± — операторы сдвига по узлам сетки на единицу вперед (+) или назад (-) по индексу т (т = г, к) от узла г, , к.

Теперь неявный оператор в левой части уравнения (46) может быть приближенно фак-торизован:

(П - С+1/2Тг - С.+-1/2Т - С+-1/2Тк 1х

X (П + С++1/2 Т+ + С-+1/2Т+ + С-+1/2Т+) Дфп+1 = ЯН5: (47)

и обращен последовательно в два дробных шага:

I шаг

= П-1(ЯН5. + С+^Дф*-. + С+ 1/2Дф*.-1,к + С+^Дф*^), (48)

II шаг

дф^1 = Дфг.д - п-1(с-1/2Дфп+1^+С-^ДФ.1^+с-^дфаи. (49)

На первом шаге по формуле бегущего счета (48) совершается разовый обход расчетной области в направлении возрастания всех индексов и определяется вспомогательная величина Дф *. На втором шаге также по формуле бегущего счета (49), но в направлении убывания всех индексов вычисляется невязка Дф:+1, по которой затем находятся значения турбулентного параметра на п + 1-м слое:

ф:+1 = ф: + Дф:+1.

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

3.1. Турбулентное течение в плоском канале

В качестве первого тестового расчета рассмотрено турбулентное течение в плоском канале (рис. 3).

Число Рейнольдса, посчитанное по половине высоты канала H и максимальной скорости течения Umax, равнялось 12300. Во входном сечении задавался равномерный профиль скорости. Длина канала принималась равной 100H. На таком расстоянии от входного сечения достигался профиль скорости, в дальнейшем не меняющийся при продвижении в продольном направлении. Расчет проводился в одной из симметричных частей канала.

Сравнение рассчитанных и экспериментально полученного в [6] профилей продольной составляющей скорости в выходном сечении канала приведено на рис. 4. На рис. 5 эти же профили приведены в координатах закона стенки и и+. Показаны параболический профиль скорости в ламинарном подслое (4) и его продолжение, а также логарифмический профиль (9) феноменологической модели турбулентного пограничного слоя. В течениях при больших числах Рейнольдса Re > 2 • 104 (по данным работы [4]) измеренные и рассчитанные профили скорости очень хорошо совпадают с логарифмическим профилем (9)

-Стенка канала Рис. 3. Расчетная область в плоском канале. П С Э1

Рис. 4. Профили продольной составляющей скорости турбулентного течения в выходном сечении плоского канала.

Рис. 5. Профиль скорости в вязком и логарифмическом слоях.

в зоне полностью развитого турбулентного пограничного слоя 30 < у+ < 400. Это означает, что область логарифмического закона является характерной чертой пристенного турбулентного пограничного слоя для указанных чисел Рейнольдса. Однако число Рейнольдса в эксперименте и в проведенном расчете равно только 12300. Поэтому наблюдается изменение закона стенки в логарифмическом слое: рассчитанные и экспериментальные значения приподняты здесь над универсальной прямой, соответствующей высоким числам Рейнольдса.

3.2. Турбулентное течение в плоском канале за обратным уступом

Следующим тестовым расчетом, на котором апробировался предложенный метод, было моделирование течения за обратным уступом. Расчетная область схематично представлена на рис. 6. Неортогональная сетка (рис. 7) состояла из 130 узлов в продольном направлении и 72 узлов в поперечном. Число Рейнольдса в данном расчете, посчитанное по входной скорости и высоте ступеньки, равнялось 37423.

Линии тока течения в окрестности уступа, характеризующие размер рециркуляционной зоны в этом месте, изображены на рис. 8. Одной из наиболее часто рассматриваемых характеристик, сравниваемых при расчетах течений данного класса, является размер от-

Рис. 6. Расчетная область течения в плоском канале за обратным уступом.

рывной зоны. В проведенных расчетах размер отрывной зоны, отнесенный к высоте уступа, равнялся 6. В экспериментально полученных по этой задаче данных [7] значение точки присоединения потока равно хг = 6.1.

Рис. 7. Фрагмент расчетной сетки в окрестности уступа.

Рис. 8. Линии тока в окрестности уступа.

Рис. 9. Распределения коэффициента давления Ср вдоль нижней стенки канала.

На рис. 9 приведено сравнение распределений коэффициента давления р вдоль нижней стенки канала, полученных численно и экспериментально.

3.3. Течение в рабочем колесе радиально-осевой гидротурбины

После проведенных на тестовых задачах сравнений численных результатов с имеющимися экспериментальными данными разработанный алгоритм был применен для расчета турбулентных течений в рабочем колесе радиально-осевой гидротурбины. Исследования данного течения в рамках невязкой модели были проведены авторами ранее и представлены в работе [1]. По своему характеру данное течение является турбулентным, и его расчет в рамках усредненных уравнений Рейнольдса с использованием некоторой модели турбулентности представляет интерес как с точки зрения уточнения качественных и количественных характеристик течений, получаемых в результате численных экспериментов, так и с методической стороны, заключающейся в возможности использования при расчетах сильно закрученных течений в элементах гидроустановок двухпараметрической к — б модели турбулентности.

На рис. 10 показаны все лопасти рабочего колеса, осуществляющего равномерное вращение вокруг своей оси с угловой скоростью На рис. 11 изображена расчетная область между двумя соседними лопастями, представляющая один из сегментов цикличности течения в рабочем колесе. Расчетная сетка имела 43 х 30 х 30 (38700) ячеек. Требовалось около 1300 шагов по времени до получения стационарного решения с точностью 10-4. Общее время решения задачи на компьютере РепИит-200 составляло порядка 1 часа.

Рис. 10. Математическая модель лопастей рабочего колеса радиально-осевой гидротурбины.

На рис. 12 представлены профили продольной составляющей скорости Cz вдоль радиуса r в сечении z = const, расположенном недалеко от выходного сечения расчетной области (см. рис. 11). Данные получены численно в рамках уравнений Эйлера [1] (пунктирная линия) и Рейнольдса (сплошная линия), а также экспериментально в лаборатории

водяных турбин Ленинградского металлического завода. Наблюдается хорошее соответствие приведенных результатов.

Рис. 11. Расчетная область между двумя соседними лопастями.

Рис. 12. Профили продольной скорости.

4. Заключение

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

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

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

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

[1] Грязин Ю. А., Черный С. Г., Шаров С. В., ШАШКИН П. А. Об одном методе численного решения трехмерных задач динамики несжимаемой жидкости. Докл. РАН, 353, №4, 1997, 478-483.

[2] DURBIN P. A. Application of a near-wall turbulence model to boundary layers and heat transfer. Int. J. Heat Fluid Flow, 14, 1993, 316-323.

[3] SHIMOMURA Y. A. Statistically derived two-equation model of turbulent shear flow in a rotating system. J. Phys. Soc. Japan, 58, 1989, 352-355.

[4] LAUNDER B. E., Spalding D. B. The numerical computation of turbulent flows. Comp. Methods in Appl. Mech. and Eng., 3, 1974, 269-289.

[5] Yakhot V., Orzag S. A., Thangam S., Gatski T. B., Speziale C. G. Development of turbulence models for shear flows by a double expansion technique. Phys. Fluids, A7, 1992, 1510.

[6] LAUFER J. Investigation of turbulent flow in a two-dimensional channel. NASA Rept., 1053, 1951.

[7] Driver D.M., Seegmiller H. L. Features of a reattaching turbulent shear layer in divergent channel flow. AIAA J., 23, No. 2, 1985, 163-171.

Поступила в редакцию ,5 октября 1998 г.