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

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

Том 6, № 5, 2001

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

__ _ и <_> , . _____

ДВУХСЛОЙНОЙ (к - е)-МОДЕЛИ*

А. В. Кузьминов, В. Н. Ллпин, С. Г. Черный Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: lapvas@ngs.ru

The method of closing for Reynolds equations, which is based on two-layer turbulent model is suggested. Two-layer model includes (k—e)-model for the fully developed turbulent flow and model with one equation for the flow near the wall. It is not nessesary to distinguish explicitly the regions, where each model is applied, because the method switches automatically from one model to another.

Введение

Новый эффективный численный метод решения трехмерных уравнений Эйлера и Навье — Стокса, записанных в простейших переменных, был предложен в [1]. В [2] он распространен на трехмерные уравнения Рейнольдса с замыкающей двухпараметрической (к — е)-моделью турбулентности. Стандартная (к — е)-модель турбулентности непригодна вблизи стенки. Известная ее модификация для низких чисел Рейнольдса (low Reynolds number (k — e)-model), решающая проблему расчета турбулентных течений в пристеночном слое посредством демпфирующих функций [3], требует сильного сгущения узлов сетки к стенке, что в трехмерном случае приводит к значительным затратам ресурсов компьютера. В большинстве практических приложений (к — е)-модели расчет турбулентных течений вблизи стенки основывается на методе пристеночных функций, использующем логарифмический закон стенки (high Reynolds number (k — e)-model). Однако при его реализации требуется, чтобы ближайший к стенке узел сетки находился в логарифмическом слое внутренней области турбулентного пограничного слоя. Вязкие эффекты должны быть здесь полностью подавлены турбулентной вязкостью. Для типичных пристеночных течений это означает, что величина безразмерного расстояния у+ от стенки до данного места должна быть больше 30. В этом случае можно использовать универсальность профиля скорости, гипотезу локального равновесия энергии турбулентных пульсаций, а также свойства локальной изотропности диссипирующих вихрей. Если же окажется, что ближайшие узлы сетки расположены в вязком или буферном подслое, т. е. 0 < У+ < 30, то метод пристеночных функций перестает здесь работать и требуется привлечение других моделей, например, алгебраической модели Прандтля.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант №01-01-00799), Программы интеграционных фундаментальных исследований СО РАН (проект №1) и Федеральной целевой программы “Интеграция” (проект 274).

© А. В. Кузьминов, В.Н. Лапин, С. Г. Черный, 2001.

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

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

1.1. Уравнения Рейнольдса

В декартовой системе координат х\, х2, х3 уравнения Рейнольдса несжимаемой жидкости с замыканием по гипотезе Буссинеска записываются в виде

дщ

дхі

0,

(1)

диі

+

д

ді дх.

_д_

дхі

_д_

дх,;

(V + Щ)

дхі дх

где щ — осредненные компоненты вектора скорости в декартовой системе координат; к — кинетическая энергия турбулентности; V — кинематическая вязкость; щ — турбулентная вязкость; р — осредненное давление.

1.2. Модели турбулентности

Для замыкания уравнений Рейнольдса (1) необходимо привлечь соотношения, позволяющие находить к и ^. При этом модель, их описывающая, должна быть относительно простой, что важно при расчете сложных пространственных течений. В то же время она

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

Наиболее подходящей с точки зрения данных требований является хорошо зарекомендовавшая себя в последнее время двухслойная модель [5, 6]. Ее основу составляют две модели: модель с одним уравнением, которая вполне адекватно описывает турбулентные характеристики в пристеночном слое, и стандартная (к—е)-модель для расчета параметров в области развитого турбулентного течения.

Модель с одним уравнением. Справедливое как в области развитого турбулентного течения, так и в вязком подслое уравнение переноса для к записывается следующим

образом [4]:

дк д ( дк\ 2vk

dt + dXj Г11 — Vk oXJ = G — e — IF' (2)

где e — скорость диссипации турбулентности; у — кратчайшее расстояние до ближайшей твердой стенки;

^ _ 0 dUi 0 _dUi | дщ _ |

G — VtSij^---, Sij — —-------+ —-, vk — v + vt-

дХу дХ1 дХг

Турбулентная вязкость определяется соотношением вида

V — CX к1/2- (3)

При вычислении e считается, что справедлива гипотеза длины пути смешения Прандтля и

_ ^ к3/2

e — CD——. (4)

1г

Входящие в (3) и (4) масштабы lv и 1г обеспечивают необходимые эффекты затухания в пристеночной области и определяются в соответствии с [7] по формулам

lv — Cl y [1 — exp ( — Rey/Av)] , 1г — Cl y [1 — exp ( — Rey/Аг)] .

Коэффициенты в формулах имеют следующие значения Ci — kar C-3/4, kar — 0.42, Аг —

2Ci, C — 0.09, Cd — 1, Av — 70, Rey — к1/2у/и.

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

I ит

У — у—,

V

где ит — sfrW есть динамическая скорость, а tw — напряжение трения на обтекаемой поверхности, должно быть больше 30. Турбулентная вязкость в этой модели определяется по формуле ,2

Vt — C, - (5)

Турбулентная кинетическая энергия находится из уравнения (2), а скорость ее диссипации — из уравнения переноса

e „ „ e2

Vt

где ve — v +---------------, аг

Ог

дe д дe '

at + дХ1 [eUj — Ve д Xj /

.3; 1. — 1.44, C& — 1.92

C£l-G — C£2T, (6)

кк

Двухслойная модель. Данный подход иногда называют “зональным моделированием”. Область течения разбивается на две зоны. Одна из моделей для высоких чисел Рейнольдса используется в основной, удаленной от стенки области течения. В пристеночной области применяется модель с одним уравнением. Это позволяет избежать недостатков (к — є)-модели с демпфирующими функциями для низких чисел Рейнольдса, одним из которых является требование иметь очень подробную сетку вблизи стенок. По точности же двухслойная модель не уступает моделям с демпфирующими функциями. Но двухслойная модель обеспечивает адекватное описание явления только в том случае, если пристеночная зона, в которой хорошо работает модель с одним уравнением, правильно выделена. Обычно граница раздела областей применимости моделей турбулентности с одним и двумя уравнениями выбирается априори, и тогда возможно переключение одной модели на другую при переходе через эту границу. Для большинства задач заранее выделить пристеночный слой затруднительно.

В настоящей работе переход от одной модели к другой осуществляется с помощью управляющей функции Л = Л (Rey), аргументом которой является турбулентное число Рейнольдса Rey. Будем считать, что на границе между вязким подслоем и областью развитого турбулентного течения имеет место значение турбулентного числа Рейнольдса Rey. Определим тогда управляющую функцию Л (Rey) следующим образом:

Л (Rey) =

1

1 + tanh

Rey — ReJ'

где А = 1 ^ 10. Число характеризует плавность перехода значений Л (Иеу) от нуля к единице в окрестности И,еу = И,е*. График функции приведен на рис. 1.

Рис. 1. График функции Л (Rey).

Уравнения для £ из модели с одним уравнением и из стандартной (к — £)-модели объединяются далее с помощью функции А в одно выражение

дє л д

dt + дх,

дє дх,

=л Í С£\—G — С£2 — j + (1 — л) í Cd і

к3/2

є

(7)

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

к2

щ = АС-------+ (1 — А) к1/2.

2

2

В пристеночной области (А = 0) уравнение переноса и диффузии для £ преобразуется в обыкновенное дифференциальное уравнение

¿£ / к3/2

Л = {С°— — £

имеющее при предположении постоянства величины Со к3/2/1£ решение

£ (¿) = Со к3/2/1£ + (£° — Со к3/2/1£) ехр(—¿).

В процессе установления решения при Ь ^ то реализуются формулы соответствующих моделей.

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

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

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

2.2. Метод решения уравнений двухслойной (к — в)-модели турбулентности

Численный метод решения уравнений стандартной (к—£)-модели подробно описан в работе

[2]. Для обеспечения консервативности и экономичности метода при дискретизации замыкающих уравнений так же, как и для основных уравнений, используется неявный метод конечных объемов. В результате получается схема, имеющая третий порядок аппроксимации членов конвективного переноса и второй — диффузионных членов. В настоящей работе уравнение для к (2) совпадает с используемым в [2] и решается таким же методом. В уравнении для £ (7) в отличие от применяемого в [2] появились дополнительные коэффициенты и источниковые члены, численное представление которых требует модификации исходного метода. Рассмотрим этот вопрос подробнее.

Перепишем уравнение (7) следующим образом:

д£

+ Ад ('£и] — у£ дх ) + ( АС£2 Т + 1 — А ) £= АС£1 т.^ + (1 — А) С0

где

Т = тах

£

к

кз/2

(8)

> 0

I

£

есть турбулентный масштаб времени. Величина є во всех членах левой части уравнения (8) при аппроксимации эквивалентного ему интегрального закона сохранения берется неявно с (п + 1)-го слоя по і. Тогда получается схема вида

єп+1 - є

^і]к

із к

л Е

т=і,з,к

Ді

дє

Уі]к + Л Е [(єи ■ Б) т+1/2 (єи ' Б)т—1/2

т=і,з,к

п+1

дє

дє

+

' дх 1

т+1/2

дє

дє

дє

»Є \ Бх^-------+ Бу^-----+ о

дх ду дг

т- 1/2

п+1

ЛСє2 Т + 1 - Л ) є^1У

із к

є к3/2\П '

ЛСє1 -о + (1 - Л) Сп—) У

+

(9)

із к

Здесь У,к — объем текущей ячейки расчетной сетки; £,к — среднее значение £ в данной ячейке; и = (и1,п2,п3) — вектор скорости; 8х,Бу,Бг — компоненты нормального к грани т +1/2 — вектора 8т+1/2, модуль которого равен площади этой грани.

Аппроксимация невязких и вязких потоков в (9) осуществляется так же, как ив [2]. Затем проводятся линеаризация и ЬИ-факторизация схемы, после чего она представляется в виде двух дробных шагов:

П 1 (уКНЗПк + С+_ 1/2Д£*_1,. + С+1/2Д£*,_ц + +-1/2^£1,к-1) >

Д,1 = Д4к — П-1 (с^Д^ + 7+1/2Д£’5+11к + ¿+1/2Д£П,+|^ .

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

О 7 Т СП

ПП 5ізк

Л

т=і,з,к

єПиП+‘-^т+1/2 - (єпип+1-^т—1/2

+

+Л

т=і,з,к

дє дє дє \

»є( Бх — + Бг

дх ду дг)

+

т+1/2

»є [ Бх дх + Бу ^ + Бг ^ дх ду дг

+

т—1/2

Л ( Сє1 к° - Сє2 Тє ) + (1 — Л) ( -є

к3/2

Ст+1/2 = 0-5Л (и ■ Б ± |и ■ Б|)т+1/2 ± Л»є

Б • Б

У

із к

т+1/2

п

Уі

із к

Ді

+ (ЛСє2 /Т + 1 - Л) Уізк + Сі—1/2 - Сі+1/2 + С7+—1/2 - С7+1/2 + С+—1/2 - С—

і—1/2 ^і+1/2 1 ^з —1/2 ^?+1/2 1 к—1/2 ^к+1/2’

п+1 п

є — є .

ДєП+1 ___ ^п+1 ¿п

п

п

п

2.3. Граничные условия

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

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

где Ск = 0.01 ^ 0.001. Скорость диссипации энергии турбулентности определяется по формуле

где l = min (0.41y, 0.095). Здесь y — расстояние до стенки, 5 — толщина пограничного слоя. Давление экстраполируется из расчетной области.

Выходная граница. На выходной границе задаются давление и поперечные составляющие скорости u2, u3, равные нулю. Для продольной компоненты скорости ui, параметров к и £ ставятся “мягкие“ условия

где £ — продольная координата расчетной сетки.

Твердые стенки. На стенках нормальная составляющая скорости полагается равной нулю. Задание нулю касательной составляющей скорости (в силу невозможности значительно уменьшать здесь шаг по нормальной к стенке координате) приводит к занижению нормального градиента продольной составляющей скорости на стенке и неправильной передаче здесь напряжения трения. Поэтому в настоящей работе для касательной составляющей скорости на стенке ставится условие частичного скольжения.

Пусть ир — вектор скорсти, рассчитанный из уравнений в ближайшем к стенке узле Р. Симметричным относительно стенки (точки т) узлу Р будет фиктивный узел Ф (рис. 2).

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

Если на стенке задавалось бы только условие непротекания (уравнения Эйлера), то необходимо было бы положить

ko = Ck U2,

dui дк д£

д£ д£ д£

u Ф = ифС = ифтС + и Фп,

/' У f'

Ф

иФтС

и фа'

иФс

Рис. 2. Граничные условия для компонент скорости на твердой стенке.

где ифтс = ирт, и.фп = — ирп. Тогда скорость на стенке иад была бы равна

иад = 0.5 (ир + иф) = ирт.

В случае задания на стенке условия прилипания (уравнения Навье — Стокса) необходимо было бы положить в фиктивном узле

иФ = иФА = иФт А + иФп,

где ифта = — ирт. Скрость на стенке была бы равна

иад = 0.5 (ир + иф) = 0.

В настоящей работе скорость в фиктивном узле задается выражением

иф = аифс + (1 — а) и фа, и скорость на стенке принимает значение

иад = 0.5 (ир + иф) = аирт.

Коэффициент а задается исходя из логарифмического закона стенки, а именно

^ . (10)

ду каг ■ ур

Напряжение трения в определении динамической скорости ит = ^/Тй полагается приближенно равным

д | ирт |

т'ш —

ду

д | и |

При подстановке значения ит = * I иь—о—~ в (10) получается следущее выражение для

11 ду

нормальной производной от касательной составляющей скорости в узле Р:

д | ирт | = у/щ д | ирт |

ду каг ■ ур у ду ’

которое преобразуется к виду

д | ирт | V

ду каг2 ■ ур ’

П й д |ирт| „ „

После аппроксимации производной —------------------------простейшей разностью первого порядка

ду

д | ирт | _ | ирт | — а | ирт |

ду ур

из (11) получается выражение для а

а = 1------V----1 Й-2 . (12)

ур | ирт | каг2

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

Кинетическая энергия турбулентности к на стенке полагается равной нулю. Для скорости диссипации ставится условие де/ду = 0.

(11)

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

Для апробации предложенного метода так же, как ив [2], рассмотрены две задачи: о течениях в плоском канале и за уступом.

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

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

График зависимости логарифма максимальной невязки Err решаемых стационарных уравнений от номера итерации N (временного слоя) приведен на рис. 4.

Рис. 3. Расчетная область течения в плоском канале.

Рис. 4. Сходимость метода на задаче о течении в плоском канале.

Сравнение профилей продольной составляющей скорости в выходном сечении канала, рассчитанных при различных значениях Ие*, а также полученного экспериментально в [9], приведено на рис. 5.

При значении Ие*, равном 100, наблюдается максимальное соответствие численно и экспериментально полученных профилей скорости. Дальнейшее увеличение Ие^ (на ри-

сунке Ие* = 250) не дает существенных изменений скорости, что свидетельствует о том, что в данной задаче модель турбулентности с одним уравнением хорошо описывает течение во всей области. Это подтверждают и рис. 6, 7.

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

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

На рис. 6 эти же профили приведены в координатах закона стенки у+ и и+. Показаны также параболический профиль скорости в ламинарном подслое и+ = у+ и логарифмический профиль и+ = -— 1п (Еу+) феноменологической модели турбулентности погра-

каг

ничного слоя. Здесь Е — постоянная, определяющая степень шероховатости стенки (для гладкой поверхности стенки экспериментально установлено Е = 8.8) [10].

На рис. 7 приведены в том же сечении профили кинетической энергии турбулентности к, обезразмеренные на и2.

Сравнение результатов расчетов осредненной величины касательной компоненты напряжения иу/и2, представленное на рис. 8, подтверждает то, что в расчете обтекания

0.1 1 10 100 юоо

У

Рис. 7. Профиль кинетической энергии турбулентности k.

2 4 г---1---1—1—1—i---1---1—1—1—i---1---1—1—1—i---1----1—1—г

uv/u2 ;

_0 2 t___I___I III____I___I III____I___I III____I___I I I

0.1 1 10 100 1000

y+

Рис. 8. Профили касательной компоненты рейнольдсовых напряжений uv/u2.

пластины двухслойная модель турбулентности с числом перехода Rey, равным 100, наиболее адекватно передает картину течения. Отметим, что уменьшение Rey до 10 в расчетах

кроме искажения параметров течения повлекло также ухудшение сходимости итераций.

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

Следующим тестовым расчетом, на котором апробировался предложенный метод, было моделирование турбулентного течения в плоском канале за обратным уступом. Расчетная область схематично представлена на рис. 9. Для рассмотрения свойств предложенного алгоритма на криволинейных сетках построена неортогональная сетка (рис. 10), состоящая из 130 узлов в продольном направлении и 72 в поперечном. Число Рейнольдса в данном расчете, посчитанное по средней входной скорости и-ш и высоте ступеньки Н, равнялось 37 423.

^ ^ / ►'Иіп > у\ /ш;п Ее = = 37 423 ' V . ' і Н <

,1 , к ^ у у у у у у / / 'оЧ

ш У У* У* 4 * 30 /г X

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

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

Линии тока течения в окрестности уступа, характеризующие область рециркуляционного течения в этом месте, изображены на рис. 11. Одной из наиболее важных характеристик, сравниваемых при расчетах течений данного класса, является размер отрывной зоны. В проведенных расчетах размер отрывной зоны, отнесенный к высоте уступа, равнялся 6, а в экспериментально полученных по этой задаче данных [11] — хг = 6.1.

На рис. 12 приведены распределения коэффициента трения

С = т'ш

7 (0.5^)

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

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

Cf 0.0020

0.0015

0.0010

0.0005

0.0000

-0.0005

-0.0010

-0.0015

-5 0 5 10 15 20 25 30 x/h

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

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

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

[2] Черный С. Г., Шашкин П. А., Грязин Ю. А. Численное моделирование пространственных турбулентных течений несжимаемой жидкости на основе (к — є)-моделей // Вычисл. технологии. 1999. Т. 4, №2. C. 74-94.

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

[4] Чжен К. Ю. Расчет течений в каналах и пограничных слоях на основе модели турбулентности, применимой при низких числах Рейнольдса / / Ракетная техника и космонавтика. 1982. Т. 20, №2. C. 30-37.

□

□ □ □

[5] JONGEN T., MARX Y. P. Design of an unconditionally stable, positive scheme for the k — £ and two-layer turbulence models // Comp. Fluids. 1997. Vol. 26, No. 5. P. 469-487.

[6] Xu D., LESCHZINER M. A., KHOO B. C., Shu C. Numerical prediction of separation and reattachment of turbulent flow in axisymmetric diffuser // Comp. Fluids. 1997. Vol. 26, No. 4. P. 417-423.

[7] WOLFSHTEIN M. The velocity and temperature distribution in one-dimentional flow with turbulence augmentation and pressure gradient // Intern. J. of Heat and Mass Transfer. 1969. Vol. 12. P. 301-318.

[8] Грязин Ю. Я., Черный С. Г., Шаров С. В. Об использовании методов типа попеременно-треугольных решения неявных разностных схем для трехмерных уравнений динамики несжимаемой жидкости // Вычисл. технологии. 1995. Т. 4, №13. C. 306-320.

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

[10] БЕЛОВ И. А., Исаев С. А., Коровков В. А. Задачи и методы расчета отрывных течений несжимаемой жидкости. Л.: Судостроение, 1989.

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

Поступила в редакцию 18 мая 2001 г., в переработанном виде — 24 июля 2001 г.