CC BY
21
УДК 532.5.011.12
В.А. Гущин1,2, П.В. Матюшин1
1 Институт автоматизации проектирования РАН 2 Московский физико-технический институт (государственный университет)
Математическое моделирование течений несжимаемой жидкости
Для прямого численного моделирования пространственных течений однородной несжимаемой вязкой жидкости используется метод расщепления по физическим факторам 8М]Р-МЕРАНЖ с явной, гибридной конечно-разностной схемой (второй порядок аппроксимации по пространственным переменным, минимальная схемная вязкость и дисперсия, монотонность). Демонстрируется использование метода 8М]Р-МЕРАНЖ для моделирования внешних и внутренних пространственных течений несжимаемой вязкой жидкости. Приводятся детальные механизмы формирования вихревых петель в следе за сферой для различных режимов течений однородной по плотности жидкости (в широком диапазоне чисел Рейнольдса 270 < Ие < 103). Рассматриваются две моды пространственного течения за «бесконечным» круговым цилиндром при 200 < Ие < 400. Демонстрируются примеры моделирования воздушных потоков в «чистой комнате».
Ключевые слова: математическое моделирование, несжимаемая вязкая жидкость, метод расщепления по физическим факторам, обтекание сферы, кругового цилиндра, воздушные потоки в «чистой комнате».
I. Введение
Многие явления, наблюдаемые в атмосфере и океане, а также проблемы гидромеханики, гидравлики, акустики, физиологии кровообращения, организации технологических процессов в силу умеренных скоростей перемещения среды можно изучать в рамках модели несжимаемой вязкой жидкости. В настоящее время практика выдвигает перед учёными-приклад-никами различного рода задачи, полное исследование которых может быть проведено в большинстве случаев лишь путём вычислительного эксперимента (ВЭ) или с помощью тщательно поставленного физического эксперимента. Однако представляющие практический интерес явления и технологические процессы либо не поддаются всестороннему физическому моделированию, либо затраты на проведение таких экспериментов чрезмерно велики.
Задачи, представляющие практический интерес, как правило, характеризуются многомерностью, нестационарно-стью, нелинейностью, наличием свободных границ и пограничных слоёв и описываются уравнениями Навье-Стокса. Нелинейность уравнений Навье-Стокса и нали-
чие малого параметра при старших производных (особенно при больших числах Рейнольдса) создают серьезные трудности как при их аналитическом исследовании (оно, по существу, возможно лишь для модельных уравнений или частных задач), так и при численном решении этих уравнений с помощью ЭВМ.
В последнее время в связи с бурным развитием вычислительной техники и вычислительной математики значительно возрасла роль ВЭ для моделирования сложных нелинейных задач гидродинамики. Основные принципы, положенные в основу ВЭ, достаточно подробно изложены в работах Н.С. Бахвалова, О.М. Белоцер-ковского, Г.И. Марчука, А.А. Самарского, П. Роуча, Н.Н. Яненко, их коллег и учеников.
Технологический цикл ВЭ включает в себя следующие основные этапы: построение физической модели — на этом этапе происходит выяснение основных физических процессов и механизмов, присущих изучаемому явлению или технологическому процессу; для получения количественной информации строится соответствующая физической математическая модель — система алгебраических, дифференциальных, интегральных или интегро-
дифференциальных уравнений с необходимыми начальными и граничными условиями; далее следует один из важных этапов в технологическом цикле ВЭ — конструирование вычислительного алгоритма, включающего два основных момента: построение дискретной математической модели, то есть аппроксимация исходной задачи разностной, и разработка эффективного метода решения разностной задачи; далее следует написание программы, набор и отладка её на ЭВМ; после этого проводится ряд тестовых расчётов, позволяющих судить о правильности и точности работы алгоритма и программы; и, наконец, находятся численные решения исходной математической задачи. В заключение проводится сравнение полученных результатов с экспериментальными данными или расчётами других авторов, на основании которого делается вывод об адекватности выбранной физической модели исследуемому явлению. При этом следует иметь в виду, что уровень используемых физической и математической моделей, алгоритм и его точность, а также конфигурация, производительность и оперативная память ЭВМ должны быть взаимно сбалансированы.
При разработке физической модели в первую очередь необходимо учесть главные особенности явления и попытаться на максимально простой модели получить адекватные результаты, отражающие суть явления.
Так, очевидно, в отрывных течениях вязкой жидкости основным является механизм, обусловленный наличием прилипания жидкости к твёрдой поверхности обтекаемого тела. Ошибки в моделировании этого механизма, имеющего молекулярную природу, могут существенно повлиять на получаемые результаты. Таким образом, при моделировании отрывных течений особое внимание следует обращать на адекватное воспроизведение течения около поверхности тела, что особенно важно с увеличением числа Рейнольдса (с уменьшением толщины пограничного слоя). Более того, в диапазоне критических чисел Рейнольдса, когда происходит потеря устойчивости ламинарного пограничного слоя и внутри пограничного слоя возникают крупномасштабные вихревые образования — пятна, необходимо правильно воспроизводить рождение и дина-
мику этих образований, то есть для расчётов таких режимов необходимо использовать конечно-разностные сетки с соответствующей разрешающей способностью. С другой стороны, течение в ближнем следе тела конечного размера как для ламинарных, так и для турбулентных режимов обтекания характеризуется наличием в нём крупномасштабных вихрей, размеры которых соизмеримы с линейными размерами тела, и, следовательно, разрешающая способность конечно-разностной сетки здесь может быть существенно ослаблена. Конечно, при этом могут быть допущены ошибки в аппроксимации диссипативного механизма таких вихрей, однако основные черты явления, включая многие количественные характеристики, воспроизводятся достаточно точно.
Большинство существующих методов решения уравнений Навье-Стокса не позволяет получать достоверные результаты при изучении свойств течений вязкой жидкости у тел сложной формы (например, при определении аэродинамических характеристик современных летательных аппаратов и подводных судов), особенно для больших чисел Рейнольдса, а также для турбулентных режимов обтекания. Достаточно точные количественные данные, сравнимые с физическим экспериментом, получены в основном для ламинарных режимов течения в двумерных стационарных задачах.
Существование же при обтекании тел таких явлений, как резкое изменение сопротивления, появление непредсказуемой по знаку боковой силы, связанных с отрывом потока и его турбулизацией, привлекают внимание как практиков, так и теоретиков. Изучение этих явлений даже на телах простой формы представляет существенный интерес. Этим объясняется большое количество экспериментальных и расчётных работ по исследованию отрывных течений около поперечно обтекаемого кругового цилиндра. В большинстве случаев течения являются существенно трёхмерными.
Математическое моделирование таких сложных течений, какими являются пространственные отрывные течения, особенно при больших числах Рейнольдса, течения со свободной поверхностью, течения стратифицированной по плотности жид-
кости, а также существенно трёхмерные течения в технологических устройствах и производственных помещениях специального назначения, предъявляет целый ряд требований к применяемым и разрабатываемым методам решения уравнений, описывающих эти течения. К таким требованиям относятся: высокий порядок аппроксимации конечно-разностных схем — второй и выше; минимальная схемная диссипация и дисперсия; работоспособность в широком диапазоне исследуемых параметров (чисел Рейнольдса и т. д.) и монотонность. Последнее свойство особенно важно при моделировании течений с областями больших градиентов гидродинамических параметров и при расчётах течений со свободной поверхностью, а также течений стратифицированной жидкости.
Разработка эффективных численных методов, обладающих перечисленными выше свойствами, и расчёт с их использованием существенно нелинейных течений несжимаемой жидкости являются весьма актуальными.
II. МЕРАНЖ — МЕтод РАсщепления по физическим
1 •• и
факторам для расчета течении Несжимаемой Жидкости
11.1. Краткий обзор
В настоящее время известно достаточно большое количество численных методов решения уравнений Навье-Стокса, описывающих течения несжимаемой вязкой жидкости. Большая часть этих методов была разработана применительно к системе уравнений, записанных относительно функции тока ф и вихря ш:
дш дф дш дф дш
"77Г + -ГГ- Т)-ТТ- -7Г- = Аф = и),
дъ ду дх дх ду
где А — оператор Лапласа, V — коэффициент кинематической вязкости. Общим недостатком этих методов является использование в том или ином виде граничного условия для вихря на твёрдой поверхности тела, которое отсутствует в физической постановке задачи. Наличие дополнительного итерационного процесса, связанного с указанным граничным услови-
ем для вихря, лимитирует скорость сходимости численных алгоритмов. Очевидно, что разностная схема, позволяющая рассчитывать течения вязкой несжимаемой жидкости без использования граничного условия для вихря на твёрдой поверхности, при всех прочих равных условиях обладает большей эффективностью.
Ограниченность методов решения (ф, ш) — системы связана также с существенным возрастанием объёма вычислений при обобщении их на случай пространственных течений, так как в этом случае, кроме трёх уравнений переноса завихренности, необходимо решать три уравнения Пуассона для составляющих векторного потенциала.
Трудности возникают и при попытках распространения указанных подходов для расчётов течений со свободной поверхностью, где необходимо ставить граничные условия для функции тока и вихря на заранее неизвестной и определяемой в процессе решения границе. Эти обстоятельства объясняют возросший в последнее время интерес к численному решению уравнений Навье-Стокса, записанных в естественных переменных ((V, р) — система или примитивные уравнения):
+ (у ■ У)и = — Ур + и А у (1)
V- V = 0
где V — вектор скорости, р — давление.
Важным этапом в развитии методов решения примитивных уравнений явился метод маркеров и ячеек — МАС, разработанный в Лос-Аламосской лаборатории [1]. Основными отличительными чертами этого метода являются: сочетание эйлерова и лагранжева подходов для описания движения жидкости, которое достигается использованием маркеров — невесомых и невзаимодействующих частиц, что существенно облегчает слежение за подвижными границами в процессе расчётов; использование разнесенного шаблона конечноразностной сетки для различных неизвестных функций. Метод МАС получил широкое распространение и к настоящему времени имеется целый ряд его модификаций.
Однако следует отметить, что в схемах типа МАС в силу выбранных разностных представлений выполнение условия прилипания приводит, с необходимостью, к
определению граничного значения вихря на твёрдой поверхности тела, удовлетворяющего условию Тома. Последнее, как известно, является условием первого порядка точности относительно шага пространственной сетки. Более того, все предложенные к настоящему времени модификации метода МАС не позволяют проводить расчёты при больших числах Рейнольдса, а конечно-разностные схемы этих подходов не обладают свойством монотонности.
П.2. Схема расщепления МЕРАНЖ
В отличие от классической схемы формального расщепления на одномерные операторы в предлагаемом подходе используется явная схема расщепления по физическим факторам [2-3]. Подобная схема расщепления используется в методе [4] и также в работе [5].
Пусть в некоторый момент времени Ьп = пт (т — шаг по времени, п — число шагов) известно поле скорости уп. Тогда схему определения неизвестных функций в момент времени ¿п+і = (п + 1)т можно представить в виде трёхэтапной схемы расщепления:
этап I:
v — Vа
-(Vа -V)Va - VVxua, (2)
этап II:
этап III:
Vn+1 _ V
-Vp,
(3)
(4)
где и = V, Б = V- V
Очевидно, что сумма уравнений I и III этапов даёт нам исходное уравнение движения (1), а уравнение (3) II этапа получается из (4) применением к последнему оператора дивергенции с учётом уравнения неразрывности (условия соленоидаль-ности поля скорости V ■ Уп+1 = 0).
Предлагается следующая физическая интерпретация приведённой схемы расщепления. На I этапе предполагается, что перенос количества движения (импульса единицы массы) осуществляется только за счёт конвенции и диффузии. Полученное таким образом поле скорости V*, вообще говоря, не удовлетворяет условию
несжимаемости (D = 0). Следует, од-
нако, отметить, что это промежуточное поле скорости имеет вполне определённый физический смысл. Действительно, если применить оператор rot к обеим частям уравнения (1), а также учесть тождество rot grad p = 0, то получим rot# = = rot vn+l = un+l, то есть уже на I этапе промежуточное поле скорости во внутренних точках исследуемой области течения даёт правильные вихревые характеристики.
На II этапе по найденному промежуточному полю скорости v с учётом условия соленоидальности вектора скорости vn+1 из решения уравнения Пуассона находится поле давления.
На III этапе предполагается, что перенос осуществляется только за счёт градиента давления (конвекция и диффузия отсутствуют).
II.3. Принцип построения конечно-разностной схемы
Прежде чем приступить в построению конечно-разностной схемы МЕРАНЖ для различных систем координат в случае двух и трёх пространственных переменных, напомним, что она должна удовлетворять следующим требованиям: высокий порядок аппроксимации — второй и выше; минимальная схемная диссипация; устойчивость в широком диапазоне чисел Рейнольдса и монотонность. Идею построения неоднородной (гибридной) схемы, удовлетворяющей перечисленным выше требованиям, рассмотрим на примере одномерного модельного линейного уравнения переноса
ft + afx = 0, a = const.
(5)
Введём равномерную по пространству сетку Пп = [хг = г ■ к, г = 0,1,2,...,п} и шаг по времени т. Определим на Пп сеточную функцию /П, совпадающую в узлах сетки с искомой функцией /. Пусть
+ а-
fi
i
а J а
г+1/2 U-1/2 ~h~
0
(6)
есть конечно-разностная аппроксимация уравнения (5). При этом рассматривается потоковый вариант метода в предположении, что а — скорость переноса некоторой субстанции /, а / — несомая субстанция.
Исследуем класс разностных схем, который можно записать в виде двухпараметрического семейства, зависящего от параметров а и в следующим образом:
СП __
Ji+1/2 —
а/п-1 +(1 - а — в)/п + в/Т+1,,
а ^ 0
а/п+2 + (1 — а — в )/п+1 + в/п, а < 0.
(7)
Выписав аналогично представление для /П—1/2 и разложив входящие в (6), (7) сеточные функции по формуле Тейлора в окрестности точки (г,п), первое дифференциальное приближение для уравнения (6) можно записать в виде
ft + afx
h то“2
-\a\(l + 2a-2{3)- —
Jx
(8)
При а = 0 и в = 0 из (6), (7) получается хорошо известная схема Годунова первого порядка аппроксимации, при а = 0, в = = 0,5 получается схема второго порядка с центральными разностями (ЦР), а при а = —0,5, в = 0 — схема второго порядка с ориентированными разностями (ОР).
Известно [6-7], что для уравнения (5) нельзя построить однородную монотонную разностную схему выше первого порядка аппроксимации. Поэтому построение монотонной схемы высокого порядка осуществляется либо на основе однородных схем второго порядка с применением операторов сглаживания, либо на основе гибридных схем с применением различных критериев переключения с одной схемы на другую в зависимости от характера решения, возможно с использованием сглаживания.
В данной работе строится гибридная (неоднородная) монотонная разностная схема, основанная на комбинации модифицированных схем с ЦР и ОР.
Исследуем схемы с ЦР и ОР с учётом сформулированных выше требований.
Рассмотрим класс модифицированных схем с ОР (МОР), то есть при в = 0. Требование минимума аппроксимационной вязкости, как легко видеть из первого дифференциального приближения (8), накладывает условие:
где 0 < С = ^ ^1, — число Куранта. Для модифицированных схем с ЦР (МЦР) (а = 0) требование минимума аппроксима-ционной вязкости накладывает условие
в — 0,5(1 - C).
(10)
Таким образом, получаемые модифицированные схемы МЦР при а = 0,
в = 0,5(1 — С) и МОР при а = —0,5(1 — С), в = 0 являются схемами второго порядка аппроксимации как по пространственной, так и по временной переменным и обладают нулевой схемной вязкостью в смысле равенства нулю коэффициента при второй производной в первом дифференциальном приближении. Очевидно, что каждая из этих схем работоспособна при любом неотрицательном коэффициенте вязкости V (числе Рейнольдса) при решении уравнения
ft + afx — vfx
(11)
Проведём теперь исследование монотонности построенных выше модифицированных схем МЦР и МОР. В соответствии с определением, введённым С.К. Годуновым, схема обладает свойством монотонности, если из того, что
AJm/2 = Ji+1 - fn > 0 для всех h следует
что fn+i/2 ^ 0 для всех i.
Записывая уравнение (6) в точке i + 1 и вычитая из полученного уравнение (6) в точке i, опуская громоздкие выкладки, для схемы МЦР имеем, что она монотонна при условии Af’+1/2 ^ С(с)AfП-l/2, где Z(C) — 0,5(1 — C)/(2 — C), а схема МОР —
при условии Afn+i/2 ^ °(C)AfП-l/2, где a(C) — 2(1 + C)/C, то есть обе схемы монотонны, если Afi+i/2 — ^АЛ-1/2, где
(12)
а
—0,5(1 — C ),
(9)
Принимая во внимание, что 0 ^ ((С) ^ 0,25, 4 ^ а (С) ^ +го, видно, что области монотонности рассматриваемых однородных схем имеют непустое пересечение. Таким образом, существует целый класс неоднородных схем, отличающихся условием переключения с одной однородной схемы на другую. Более общий анализ с учётом знаков коэффициента а, первых /П+1/2 и вторых разностей
А217+1/2 = А/п+1 — А/7 ПозвоЛЯет получить следующее условие переключения: если
(аА/А2/>7+1/2 > 0,
то используется схема МОР, то есть в (6):
-П
г-1 ,
rn ________
J i+1/2 =
0,5(3 - C)fn - 0,5(1 - C)fn
a ^ 0;
0,5(3 - C)f+i - 0,5(1 - C)fi+2,
a < 0.
(13)
для fn-1/2 — индекс i в (13) уменьшается на единицу, а при
(aAf Д2f)n+i/2 < 0,
используется схема МЦР, то есть
t = 15
III. Некоторые примеры рассматриваемых задач
rn ____
J i+1/2 =
= 0,5(1 - C sign a)fn+i +0,5(1 + Csign a)fn.
Построенная таким образом гибридная конечно-разностная схема (6) для уравнения (5) удовлетворяет перечисленным выше требованиям: на гладких решениях имеет второй порядок аппроксимации по временной и пространственной переменным; обладает минимальной схемной вязкостью и диссипацией; устойчива при выполнении условия Куранта 0 < C ^ 1 и монотонна.
В качестве тестового примера, демонстрирующего качества построенной гибридной схемы, рассмотрим для уравнения (5) решение задачи Коши с начальными условиями
f (0 х) = / 1 если x ^ 0 J ( , ) | 0, если x > 0;
при а = 1, h = 0,2, C = 0,5.
Для t = 15 (n = 150) на рис. 1 показана функция f П, полученная по схемам первого [6] — кривая 1, второго [8] — кривая 2, третьего [7] — кривая 3 порядков точности и с помощью предложенной гибридной схемы — кривая 4.
Отметим, что монотонная схема первого порядка точности обладает слишком большой диссипацией, в то время как предлагаемая схема, будучи монотонной, по диссипативным свойствам сравнима со схемами второго и третьего порядков точности.
III.1. Обтекание сферы
Постановка задачи. Рассматривается задача об обтекании сферы равномерным на бесконечности потоком однородной несжимаемой вязкой жидкости, движущимся со скоростью W0. Система уравнений Навье-Стокса, описывающих течение, имеет вид
dv _ 2
— + {v • V)v = —Vp + ^-Av, dt Re
V-v = 0, (14)
где v — это нормированный (на Wo) вектор скорости, p — нормированное (на pW02) давление, р — плотность жидкости.
Пусть х, y, z — декартовы координаты и (U, V, W) — составляющие вектора скорости v вдоль направлений х, y, z соответственно. Предполагается, что вектор скорости набегающего потока (0, 0, W0) параллелен оси z. Система уравнений Навье-Стокса (14) записывается в сферической системе координат R, 0, у (х = R sin 0 cos у, y = R sin 0 sin у, z = R cos 0). Для корректного моделирования отрывных течений необходимо сгущать конечно-разностную сетку к поверхности тела (разрешать пограничный слой). Поэтому в радиальном направлении вводится следующее преобразование:
R¿ = R {гі) = 1 + -Tj-N0
i = 1 : N
(15)
где No — число ячеек расчётной сетки, помещённых в пограничный слой
5 = л/2/Re. Для всех расчётов настоящей работы: m = 3, rmax = 3, RN æ 14,5d. Таким образом, исследуемая расчётная область отображается на прямоугольный параллелепипед [(г,0,ф): 0 ^ r ^ rmax, 0 ^ 0 ^ п, 0 ^ ф ^ 2п] с равномерной сеткой (N х M х L) внутри него.
Пусть и, v, w — составляющие вектора скорости v вдоль направлений R, 0, ф, соответственно. На твёрдой поверхности сферы ставятся граничные условия прилипания v = 0, w = 0 и непротекания и = 0. На внешней границе расчётной области (сферической поверхности с радиусом Rn (15)) при z < 0 ставятся граничные условия невозмущённого потока: и = cos 0, v = — sin 0, w = 0; а при z > 0: и = cos 0, v = — sin 0, dw/dR = 0.
Для моделирования отрывных течений однородной несжимаемой вязкой жидкости, описываемой уравнениями (14), используется метод расщепления по физическим факторам SMIF-МЕРАНЖ с явной, гибридной конечно-разностной схемой для аппроксимации конвективных членов уравнений, обеспечивающей второй порядок аппроксимации по пространственным переменным, минимальные схемные вязкость и дисперсию, работоспособность в широком диапазоне чисел Рейнольдса, монотонность.
Эффективность метода моделирования SMIF-МЕРАНЖ и возросшая мощность суперкомпьютеров позволила адекватно моделировать пространственные отрывные течения несжимаемой вязкой жидкости около сферы при умеренных числах Рейнольдса. Вычисления проводились на многопроцессорном вычислительном комплексе МВС-1000 (42 процессора Intel Xeon, 1.7 GHz).
Визуализация вихревых структур в следе за сферой. Для понимания динамики и механизмов отрыва вихревых структур несжимаемой вязкой жидкости в следе за сферой не достаточно трёх общеизвестных интуитивных индикаторов вихрей (минимум давления, линии тока и изоповерхности завихренности). Например, для простейшего случая осесимметричного обтекания сферы при Re = 100 в следе
можно выделить две вихревые структуры: кольцо в рециркуляционной области и оболочку вокруг кольца [9-12]. В то же время линии тока в системе отсчёта, связанной со сферой, визуализируют только кольцо, изолинии завихренности ш = 0,5го^ -только оболочку, а изолинии давления не визуализируют ни того, ни другого [9-12]. Аналогично, для сложного пространственного течения изоповерхности модуля завихренности показывают не все вихревые структуры течения, а только часть из них (рис. 2а). Поэтому были предложены специальные методы визуализации, позволяющие идентифицировать большую часть вихрей следа. Опишем вкратце суть двух методов [13, 14], которые используются в этой работе (рис. 2б, в).
а)
г) д)
Рис. 2. Двухнитевой след, Re = 250: а—в — изоповерхности |w| = 0,2(w = 0,5 • rot v), Л2 = —2 • 10_5 и в = 0,005; г-д — изолинии в > 0 (шаг 0,06) и Л2 < 0 (шаг 0,02) на фоне линий тока в плоскости симметрии следа (сетка 180 х 90 х 180)
Зафиксируем некоторую точку в поле течения и исследуем, как ведут себя линии тока в декартовой системе координат x = (x,y,z) с началом в этой точке и движущейся со скоростью этой точ-
ки. В окрестности этой точки (0,0,0) в линейном приближении можно записать V = дх/д.Ъ ~ Ох, где G — тензор градиента скорости (О^ = Vij = дvi/дxj). Если два собственных значения а1 и а2 тензора О — комплексно-сопряженные, то есть а1 = а—гв, а2 = а+гв (в = ^(а1>2) > 0), то из курса обыкновенных дифференциальных уравнений известно, что могут быть выбраны соответствующие им два комплексно-сопряженных собственных вектора 0,5(Ьх ± гЬ2), составленные из двух действительных векторов и Ь2, образующих плоскость, фазовыми траекториями в которой будут либо замкнутые овальные траектории с центром в зафиксированной точке (при а = 0), либо спирали с фокусом в этой точке, что означает наличие вихря в этой точке (см. [15], с. 119-122). При этом угловая скорость вращения жидкости вокруг этой зафиксированной точки равна в. Поэтому было введено определение сердцевины вихревого течения как совокупности подобластей течения с комплексно-сопряженными собственными значениями тензора градиента скорости [13] (см. изолинии в > 0 на рис. 2г и изоповерхность в = 0,005 на рис. 2в, где в силу плоскостной симметрии течения показана только одна её половина).
В статье [14] вводится другое определение сердцевины вихревых течений для несжимаемой жидкости и показывается его работоспособность для визуализации некоторых реальных пространственных течений несжимаемой жидкости. В [14] сердцевина вихревого течения определяется как совокупность подобластей течения с отрицательным вторым собственным значением тензора Б2 + П2 в них (Л2 < 0, рис. 2б, д), где Б и
П — это симметричная и антисимметричная части тензора градиента скорости (Б^ = 0,5(^ + Vji) — тензор скоростей деформаций; Пу = 0,5(^ , j — Vj, 4) — тензор завихренности; у симметричного тензора Б2 + П2 — три действительных собственных значения А1 ^ Л2 ^ Аз).
Для ПЛОСКИХ течений (71,2 = по"
этому определения сердцевины вихревого течения, предложенные в [13, 14], эквивалентны. В то же время для осесимметричных течений (в следе за сферой) изоповерхность в = 0 охватывает больший объём, чем изоповерхность Л2 = 0, но топо-
логии этих изоповерхностей близки. При этом экстремумы функций А2 и в в рециркуляционной области примерно совпадают с особой точкой (центром) в картине линий тока в системе отсчёта, связанной со сферой.
Режимы течений в следе за сферой.
а) Прямолинейный двухнитевой след (200 < Ив ^ 270). При
20,5 < И,е ^ 200 в следе за сферой формируется осесимметричная стационарная рециркуляционная область (в картине линий тока в системе отсчёта, связанной со сферой), которой соответствует вихревое кольцо в картинах изолиний в > 0 и Л2 < 0. В свою очередь с вихревой оболочкой следа в этих же картинах соотносится рециркуляционная зона в картине линий тока в системе отсчёта, связанной с набегающим потоком. При возрастании И,е от 20,5 до 200 кольцо укрупняется; вихревая оболочка приближается к кольцу, вытягивается вдоль линии движения сферы и соединяется с ним [10, 11]. При этом увеличивается скорость частиц жидкости в кольце и формируется локальный минимум давления в рециркуляционной области. При И,е = 200 локальные минимум давления, максимум в и минимум Л2 в рециркуляционной области примерно совпадают между собой и с особой точкой (центром) в картине линий тока в системе отсчёта, связанной со сферой. После введения в рассчитанное осесимметричное течение малого кратковременного возмущения в виде сдвига набегающего на сферу потока вдоль вертикальной оси х (дШ/дх > 0) при И,е ^ 200 течение остаётся осесимметричным, а при И,е > 200 давление вдоль сердцевины вихревого кольца перераспределяется таким образом, что в нижней его части оно становится меньше, чем в верхней. При этом молекулы жидкости в сердцевине кольца движутся по спирали сверху-вниз, а на его периферии — снизу-вверх.
Далее кольцо наклоняется и устойчивый верхний фокус на рис. 3а (соответствующий верхней части кольца) начинает подпитываться жидкостью из набегающего потока, прошедшей под сферой. (В устойчивом фокусе линия тока скручивается к его центру, а в неустойчивом — раскручивается от центра фокуса.) Из верхнего фокуса на рис. 3а жидкость идёт далее
по сердцевине кольца в нижний неустойчивый фокус, потом, раскручиваясь, уходит вниз по течению, огибая верхний фокус слева и сверху. Линии тока, отличные от линии тока на рис. 3а, берущие своё начало из набегающего на сферу потока и проходящие около линии первичного отрыва (на сфере), огибают сердцевину кольца по спирали несколько раз, поднимаясь вверх, и уходят вниз по течению, образуя две вихревые нити (рис. 2б, в). Следовательно, можно утверждать, что эти нити связаны с кольцом.
Рис. 3. Двухнитевой след, Ие = 250: а — пространственная линия тока, стартующая в плоскости симметрии следа; б — линии трения на тыльной части поверхности сферы; в — изолинии давления (шаг 0,01) на фоне линий тока в плоскости симметрии следа
Рис. 4. Цепочка вихревых петель в следе за сферой при Ие = 280; изоповерхности в = 0,02 в течение периода: а—б — Ь = 1496,4, 1504,4 (сетка 120 х 60 х 120)
При 200 < Ие ^ 211 на стадии наклона вихревого кольца процесс трансформации следа прекращается. При 211 < И,е ^ 270 этот процесс продолжается (происходят затухающие азимутальные колебания в рециркуляционной области). Более детально о характере этих колебаний будет сказано ниже. Таким образом, при 200 < И,е ^ 270 осесимметричный след трансформируется в стационарный прямолинейный двухните-вой след (рис. 2, 3), характеризующийся
наличием плоскости симметрии и отличными от нуля коэффициентами суммарной боковой силы Сг и вращательных моментов.
б) Нестационарный «волнообразный» двухнитевой след
(270 < Ив ^ 290). Как уже было отмечено выше, при 200 < И,е ^ 270 через некоторое время после введения возмущения в рассчитанное осесимметричное течение устанавливается равновесие между количествами жидкости, подпитывающей верхний фокус и проходящей по сердцевине деформированного кольца и далее вниз по течению (рис. 3а). При И,е > 270 такого равновесия не наступает, а наблюдается установившийся периодический процесс формирования новых звеньев цепочки вихревых петель в следе за сферой (рис. 4, 5а). Происходит периодическое вытягивание фиксированного (верхнего) края вихревой оболочки и отрыв края этой оболочки, индуцированное вихревыми нитями. В рециркуляционной зоне периодически генерируются новые вихревые кольца (или полукольца) и соответствующие им новые вихревые нити, которые индуцируют вытягивание фиксированного края вихревой оболочки. Так, в общих чертах можно охарактеризовать периодический процесс формирования новых звеньев цепочки вихревых петель в следе за сферой. В наших работах [16-17] этот процесс впервые был рассмотрен во всех деталях. Были выделены основные механизмы формирования вихрей, которые работают в трёх различных областях течения (рециркуляционной области, вихревой оболочке и внешнем течении). Процесс формирования новых звеньев цепочки вихревых петель был закодирован через последовательность этих основных механизмов формирования вихрей. Было показано, что при 270 < И,е ^ 290, 290 < И,е ^ 320 и 320 < И,е ^ 400 реализуются различные механизмы формирования вихрей. Различия наблюдаются в рециркуляционной области. Таким образом, при 290 < И,е ^ 320 обнаружен ранее неизвестный режим течения [16].
При И,е > 400 происходит периодическое вытягивание противоположных краев вихревой оболочки, сопровождаемое вращением этой оболочки вокруг линии движения сферы (рис. 5б). Характеристики
рассчитанных течений, такие, как усреднённые по времени коэффициенты сопротивления {Са) и суммарной боковой силы {Сг), углы отрыва, период отрыва хорошо согласуются с экспериментальными данными и расчётами других авторов. Ниже подробно рассматриваются процессы формирования новых звеньев цепочки вихревых петель в следе за сферой в течение одного периода при И,е = 280, 320, 350 и Ие > 400.
которого при Ь = 1498,4 исчезает, а при Ь = 1498,9 появляется вновь (рис. 6, 7б).
Рис. 5. Цепочки вихревых петель в следе за сферой: а) Ие = 350 (изоповерхность Л2 = —2 • 10—5); б) Ие = 700 (Л2 = —10—5)
Начнём с Ие = 280 ({Са) = 0,683, {Сг) = 0,079) и с момента времени
Ь = 1496. 4, при котором Сг близок к своему минимальному значению. На рис. 4а в вихревой структуре течения можно выделить оболочку (окружающую рециркуляционную область), две основные петли (ориентированные вверх) и одну индуцированную петлю (направленную вниз). Надо проследить, как в течение одного периода появляются новые индуцированная и основная петли. Сначала две ножки первой основной петли (две нити, выходящие из рециркуляционной области и связанные с её кольцом) генерируют во внешнем течении головную часть новой индуцированной петли и вытягивают вниз по течению боковые края оболочки (рис. 6б (3Ь и 2э)). В то же время в рециркуляционной области из-за сдвиговой неустойчивости Кельвина-Гельмгольца зарождается новое кольцо И1 (рис. 6, 7а), верхняя часть
Рис. 6. Ближний след при Ие = 280; изоповерхности в = 0,005 в течение периода: а—г — Ь = 1496,4, 1500,4, 1504,4, 1507,4; помечены области течения, где работают основные механизмы формирования вихрей: Н, 2э, 21, 31/Ь
Далее обновленное новое кольцо развивается и встает на место постепенно исчезающего предыдущего кольца рециркуляционной области (рис. 6, 7в-г). Этот процесс в картине мгновенных линий тока характеризуется увеличением верхнего фокуса (верхней части нового кольца) за счёт дополнительной подпитки его жидкостью,
прошедшей над сферой (рис. 7г, а). В то же время развитие обновленного нового кольца сопровождается образованием новых нитей, связанных с ним (рис. 6в (Н)), и генерацией во внешнем течении (над ножками первой петли) зародыша головной части нарождающейся (новой основной) петли (рис. 6в (31)). По мере развития новых нитей, соединённых с новым кольцом, связь ножек первой петли с исчезающим предыдущим кольцом рециркуляционной области пропадает (рис. 6в). Затем верхний край оболочки свертывается, отрывается и исчезает (рис. 6г (21), 7в-г), и завершается формирование новой основной петли (рис. 6г, а).
В описанном периодическом процессе можно выделить пять основных механизмов формирования вихрей, которые работают в трёх различных областях течения (рециркуляционной области ф1), оболочке ф2) и внешнем течении ф3)):
1к) генерация кольца (или полукольца) в D1 у поверхности сферы из-за сдвиговой неустойчивости Кельвина-Гельмгольца (рис. 6, 7а);
Н) образование в D1 двух нитей (ножек нарождающейся петли), связанных с этим кольцом (рис. 6в-г);
2э) вытягивание боковых краёв оболочки D2 вниз по течению (рис. 6б);
21/Ь) свертывание и отрыв верхнего или нижнего краёв D2 (рис. 6, 7в-г);
31/Ь) генерация в D3 головной части петли, ориентированной вверх или вниз (рис. 6в, б);
(процессы 2э, 21, 2Ь, 31, 3Ь индуцированы нитями, выходящими из D1).
Таким образом, детальный механизм формирования вихрей в следе за сферой в течение одного периода при 270 < Ие ^ 290 можно символически записать как М{1к-2э-1к-3Ь}-{Н-31-21} (где фигурными скобками выделены механизмы, происходящие в первой и второй половинах периодов).
Вернемся теперь к прямолинейному двухнитевому следу, а точнее к процессу его формирования. Если при 200 < Ие ^ 211 после введения возмущения осесимметричное кольцо деформируется, и структура течения далее не меняется со временем, то при 211 < И,е ^ 270 после первоначальной деформации осесимметричного кольца в рециркуляцион-
ной области генерируется ряд новых колец. Этот периодический процесс затухает со временем и его можно описать как М{1к-2в-1к-3Ь}-{Н-31} (так как не наблюдается отрыв верхнего края оболочки). Продолжительность этого процесса установления растёт с увеличением И,е (при И,е = 250 зарождается три новых кольца, а при И,е = 270 — 35 колец).
г)
Рис. 7. Мгновенные линии тока и изолинии в > 0 (шаг 0,04) при Ие = 280 в плоскости симметрии следа: а—г — Ь = 1496,4, 1500,4, 1504,4, 1507,4
в) Детальный механизм формирования вихрей при Ив = 350.
Рассмотрим теперь детальный механизм формирования новых индуцированной и основной петель в течение одного периода при И,е = 350 (рис. 5а, 8, {Са) = 0,627,
{Сг) = 0,068). Здесь также начнем рассмотрение с такого момента времени Ь = 882, при котором Сг близок к своему минимальному значению. В силу схожести на большей части периода механизмов формирования петель при Ие = 320 и 350 используем рис. 9 (Ие = 320) для дополнительной иллюстрации процессов, происходящих при Ие = 350.
Рис. 8. Ближний след при Ие = 350; изоповерхности в = 0,02 в течение периода: а-г — Ь = 882, 886,5, 889,5, 891 (сетка 180 х 90 х 180); помечены основные механизмы: Ы, И, 2я, 21, 31
Л О
.. в) г)
Рис. 9. Ближний след за сферой при Ие = 320; изолинии в > 0 (шаг 0,04) в плоскости
симметрии следа в течение периода: а-г — Ь = 453, 455, 456,5, 460 (сетка 120 х 60 х 120)
При Ь = 882 зарождается новое кольцо И1 (рис. 8, 9а (1к)), верхняя часть которого исчезает уже при Ь = 883,5. При Ие = 350 размеры верхнего и нижнего фокусов в картине мгновенных линий тока в плоскости симметрии следа естественно увеличиваются по сравнению с Ие = 280 (рис. 7а). Кроме того, при Ь = 882 нет притока в рециркуляционную область жидкости, прошедшей под сферой, что связано с формированием в этой области второго нового кольца И2 (рис. 9б), соединённого с ножками первой основной петли. Эти ножки формируют новую индуцированную петлю, вытягивая боковые края оболочки вниз по течению (рис. 8а (2в)). При этом ножки новой индуцированной петли связаны с основным кольцом рециркуляционной области (рис. 8б).
Далее верхняя часть кольца И2 генерирует третье новое полукольцо И3 в верхней части рециркуляционной области и исчезает (рис. 8б, 9в). Затем полукольцо И3 соединяется с оболочкой (рис. 8в), а ножки индуцированной петли формируют во внешнем течении зародыш головной части нарождающейся (новой основной) петли и вытягивают боковые края оболочки вниз по течению (рис. 8в-г (31 и 2в)); а верхняя часть основного кольца рециркуляционной области расщепляется на две части, левая часть сдвигается к сфере и соединяется с оболочкой (рис. 8в-г (1^), превращаясь в четвертое новое кольцо. При этом ножки новой индуцированной петли связаны уже с правой частью основного кольца, а разви-
тие четвёртого нового кольца сопровождается образованием нитей (ножек нарождающейся (новой основной) петли), связанных с четвертым кольцом (рис. 8г (Н)). В то же время наблюдается свертывание и отрыв верхнего края оболочки (рис. 8г (21)).
Таким образом, детальный механизм формирования вихрей в следе за сферой в течение одного периода при Ие = 350 можно символически записать как
М {{1к) — 2в — 1к — {1к)} —
— {{2в) — Ы — {1/ — 3Ь — 2Ь)}
(в угловых скобках помещены те механизмы, которые наблюдались и при Ие = 280). Здесь появился новый (шестой) основной механизм формирования вихрей:
Ы) сдвиг основного кольца рециркуляционной области к сфере (рис. 8в-г).
Несмотря на многочисленные сложные процессы, происходящие при Ие = 350 в рециркуляционной области, картина мгновенных линий трения на тыльной части поверхности сферы слабо меняется в течение периода. Происходят лишь незначительные колебания задней критической точки (рис. 3б). Проведённые методические расчёты при Ие = 350 на различных сетках 60 х 36 х 72, 120 х 60 х 120 и 180 х 90 х 180 показали сохранение топологии течения и детального механизма формирования вихрей.
г) Три нестационарных режима течений при 270 < Ив ^ 400. При увеличении числа Рейнольдса с 270 до 400 механизм формирования вихрей сильно изменяется. При Ие = 290 в зачаточном состоянии в течение ДЬ & 1 наблюдает-
ся существование второго нового кольца. При дальнейшем возрастании числа Рейнольдса (вплоть до Ие = 320) верхняя часть второго нового кольца увеличивается и генерирует третье новое кольцо, которое превращается потом в основное кольцо (рис. 9в-г). При этом механизм 3Ь заменяется на 2э, так как индуцированная петля при 270 < Ие ^ 290 начинает формироваться со своей головной части, а при 290 < Ие ^ 320 — с ножек. При Ие = 320 в зачаточной стадии уже наблюдается сдвиг левой части основного кольца к сфере (механизм Ы) (рис. 9г).
Таким образом, механизмы формирования вихрей при 270 < Ие ^ 290,
290 < Ие ^ 320 и 320 < Ие ^ 400 существенно различаются. Это различие наблюдается исключительно в рециркуляционной области. Процессы формирования вихрей в рециркуляционной области с большим трудом поддаются визуализации в экспериментах. Поэтому в работах [18-20] при 270 < Ие ^ 400 выделяются только два режима течения (при 270 < Ие ^ 290 и при 290 < Ие ^ 400). Иными словами, при 290 < Ие ^ 320 обнаружен новый режим течения, который характеризуется периодическим формированием второго нового кольца в рециркуляционной области. Это кольцо связано с ножками первой основной петли следа, а ножки формирующейся индуцированной петли связаны с основным кольцом рециркуляционной области. Следовательно, начиная с Ие = 290, периодически наблюдаются четыре нити, выходящие из рециркуляционной области (рис. 8б); тогда как при 270 < Ие ^ 290 визуализируются только две нити (рис. 6). При 320 < Ие ^ 400 периодически происходит сдвиг части основного кольца к сфере; а при дальнейшем увеличении числа Рейнольдса индуцированные петли становятся «двойниками» основных петель (механизмы их формирования становятся одинаковыми, то есть в течение периода реализуется механизм
М {{2в) — Ы — 1к — {1/ — 3Ь — 2Ь)} —
— {{2в) — Ы — 1к — {1/ — 3Ь — 2Ь)}).
111.2. 2Ю-3Ю переход при обтекании кругового цилиндра
Рассматривается задача об обтекании кругового цилиндра конечной длины Ь (с периодическими граничными условиями на его торцах) равномерным на бесконечности потоком однородной несжимаемой вязкой жидкости. Таким образом, симулируется обтекание «бесконечного» кругового цилиндра. Для этой задачи впервые были найдены максимальные значения разности фаз вдоль оси цилиндра, приблизительно равные 0,1Т и 0,02Т для двух мод течения А (200 < Ие < 300, 3,5й ^ А ^ 4^, рис. 10) и В (300 < Ие < 400,
0,8^ ^ А ^ 1,0й, рис. 11) соответственно [9, 12], где Т — период, А — поперечный (вдоль оси цилиндра) размер периодических вихревых структур следа.
Рис. 10. Две изоповерхности потоковой компоненты завихренности для Ие = 230, Ь = 7,5й (Мода А)
О
Рис. 11. Две изоповерхности потоковой компоненты завихренности для Ие = 320, Ь = 7,5й (Мода В)
ІІІ.3. Пример расчёта воздухопереноса в «чистой комнате»
Рассматривается задача, направленная на решение фундаментальных проблем разработки эффективных методов и средств проектирования, контроля и управления воздушными и тепловыми потоками в «чистых комнатах», предназначенных для производства изделий микроэлектроники, фармацевтики, медицины, биотехнологии и используемых в таких областях промышленности, как космическая, химическая, пищевая, и т. п, а также в ожоговых центрах и родильных домах. Используемая для математического моделирования система уравнений приведена в [17].
Пример моделирования воздушных потоков в действующей «чистой комнате»
«Озон-1», предназначенной для производства изделий микроэлектроники, приведён на рис. 12. Скорость вертикальной подачи воздуха через всю поверхность потолка составляет 0,45 м/с. Скорости откачки воздуха из одиннадцати отверстий под фальшь-полом различны (рис. 12а). Математическое моделирование показало, что отклонение линий тока от вертикали в «чистой комнате» «Озон-1» превышает предельно допустимые значения (рис. 12Ь). После проведения нескольких численных экспериментов для различных расположений отверстий в полу, было показано, что при помощи перемещения двух отверстий (рис. 12с) отклонение линий тока от вертикали во всей «чистой комнате» «Озон-1» становится удовлетворительным (рис. 12й). Таким образом, при помощи минимальных изменений конструкции пола «чистой комнаты» «Озон-1», удалось существенно улучшить её характеристики.
FIRST CONFIGURATION
С) SECOND
Рис. 12. «Чистая комната» «Озон-1»: а) первоначальная конфигурация, Ь) линии тока в вертикальной плоскости для первоначальной конфигурации комнаты; с) изменение конфигурации расположения дырок в полу; ё) линии тока в вертикальной плоскости для второй конфигурации расположения дырок
Пример моделирования воздушных потоков в предварительном проекте «чистой комнаты» «Озон-2», предназначенной для производства изделий микроэлектроники, приведён на рис. 13. Кроме сильного наклона линий тока в вертикальной плоскости (рис. 13Ь), отмеченной красной лини-
ей на рис. 13а, здесь наблюдается сильный вихрь около левой стены комнаты (выдувание пыли из под пола, рис. 13с-^. Это говорит о непригодности данного проекта «чистой комнаты» для реализации.
a) ROOM CONFIGURATION
Рис. 13. «Чистая комната» «Озон-2»: а) проектная конфигурация, Ь) линии тока в вертикальной плоскости; с-ё) линии тока в вертикальной плоскости около левой стены комнаты (выдувание пыли из под пола)
IV. Заключение
Подробно описан метод расщепления по физическим факторам БМШ-МЕРАНЖ с явной, гибридной конечно-разностной схемой (второй порядок аппроксимации по пространственным переменным, минимальная схемная вязкость и дисперсия, монотонность), который используется в настоящей работе для прямого численного моделирования внешних и внутренних пространственных течений несжимаемой вязкой жидкости. Рассмотрено несколько различных задач.
В задаче об обтекании сферы равномерным на бесконечности потоком однородной несжимаемой вязкой жидкости
приводятся детальные механизмы формирования вихревых петель в следе за сферой для различных режимов течений (в широком диапазоне чисел Рейнольдса 270 < Re < 103).
В задаче об обтекании «бесконечного» кругового цилиндра равномерным на бесконечности потоком однородной несжимаемой вязкой жидкости демонстрируются две моды пространственного течения при 200 < Re < 400.
В задаче моделирования воздухо-, тепло- и массопереноса в «чистых комнатах» продемонстрирован пример успешного использования математического моделирования для кардинального улучшения характеристик действующей «чистой комнаты» (ЧК) «Озон-1», предназначенной для производства изделий микроэлектроники. Также показан пример тестирования работоспособности ЧК «Озон-2» на стадии проекта.
Работа выполнена при поддержке Российского фонда фундаментальных исследований (гранты 08-01-00662, 08-01-91306, 08-07-00159 и 09-01-92102) и программы № 3 фундаментальных исследований ОМН РАН.
Литература
1. Harlow F.H., Welch J.E. Numerical
calculation of time: dépendent viscous
incompressible flow of fluid with free surface // Phys. Fluids. — 1965. — V. 8, N. 12. — P. 2182-2189.
2. Белоцерковский О.М., Гущин В.А., Щенников В.В. Метод расщепления в применении к решению задач динамики вязкой несжимаемой жидкости // Ж. вы-числ. матем. и матем. физ. — 1975. — Т. 15, № 1. — С. 197-207.
3. Гущин В.А. Пространственное обтекание трёхмерных тел потоком вязкой жидкости // Ж. вычисл. матем. и ма-тем. физ. — 1976. — Т. 16, № 2.
С 529-534.
4. Amsden А.А., Harlow F. Н. The SMAC method / / Los Alamos Scient. Lab. Rept. — 1970. — LA-4370.
5. Fortin M., Peyret R., Temam R. Resolution numeriques des equations de Navier-Stokes pour un fluide incompressible // J. de Mecanique. — 1971. — V. 10, N. 3. — P. 357-390.
6. Годунов С.К. Разностный метод численного расчёта разрывных решений гидродинамики // Матем. сб. — 1959.
Т. 47. — С. 271-306.
7. Холодов А.С. О построении разностных схем с положительной аппроксимацией для уравнений гиперболического типа // Ж. вычисл. матем. и матем. физ. — 1978. — Т. 18, № 6. — С. 1476-1492.
8. MacCormack R.W. The effect of viscousity in hypervelocity impact cratering // AIAA Paper. — 1969. — 2, 69-354. -P. 1-7.
9. Gushchin V.A., Kostomarov A.V., Matyushin P.V., Pavlyukova E.R. Direct numerical simulation of the transitional separated fluid flows around a sphere and a circular cylinder // Jnl. of Wind Engineering
& Industrial Aerodynamics. — 2002.
V. 90/4-5. — P. 341-358.
10. Гущин В.А., Матюшин П.В. Классификация режимов отрывных течений жидкости около сферы при умеренных числах Рейнольдса // Математическое моделирование: проблемы и результаты. — М.: Наука, 2003. — С. 199-235.
11. Матюшин П.В. Численное моделирование пространственных отрывных течений однородной несжимаемой вязкой жидкости около сферы: Дис. на соиск. степени канд. физико-математических наук: 05.13.18. — М., 2003. — 194 с.
12. Gushchin V.A., Kostomarov A.V., Matyushin P.V. 3D Visualization of the Separated Fluid Flows / / Journal of
Visualization. — 2004. — V. 7, N. 2. — P. 143-150.
13. Chong M. S., Perry A.E., Cantwell B.J. A general classification of three-dimentional flow field // Phys. Fluids. — 1990. — V. A2. — P. 765-777.
14. Jeong J., Hussain F. On the identification of a vortex // J. Fluid Mech. — 1995. — V. 285. — P. 69-94.
15. Понтрягин Л.С. Обыкновенные дифференциальные уравнения. — М.: Наука, 1974. — 331 c.
16. Гущин В.А., Матюшин П.В. Ме-
ханизмы формирования вихрей в следе за сферой при 200 < Re < 380//
Изв. РАН. Механика жидкости и газа. — 2006. — № 5. — С. 135-151.
17. Гущин В.А., Матюшин П.В. Математическое моделирование пространственных течений несжимаемой жидкости // Математическое моделирование. — 2006. — Т. 18, № 5. — С. 5-20.
18. Magarvey R.H., Bishop R.L. Transition ranges for three-dimensional wakes // Can. J. Phys. — 1961. — V. 39. — P. 1418-1422.
19. Sakamoto H., Haniu H. A study on vortex shedding from spheres in a uniform flow // Trans. ASME: J. Fluids Engng. — 1990. — V. 112. — P. 386-392.
20. Sakamoto H., Haniu H. The formation mechanism and shedding frequency of vortices from a sphere in uniform shear flow // J. Fluid Mech. — 1995. — V. 287. — P. 151-171.
Поступила в редакцию 15.09.2009.
