РАСЧЕТ ДВИЖЕНИЯ ВЯЗКОЙ ЖИДКОСТИ, ЧАСТИЧНО ЗАПОЛНЯЮЩЕЙ ВРАЩАЮЩУЮСЯ ПОЛОСТЬ, ПРИ БОЛЬШИХ ЧИСЛАХ РЕЙНОЛЬДСА Текст научной статьи по специальности «Физика»

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

Том 24, № 3, 2019

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

Ю. В. ПИВОВАРОВ

Институт гидродинамики им. М.А. Лаврентьева СО РАН, Новосибирск, Россия Контактный e-mail: pivov@hydro.nsc.ru

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

Ключевые слова: конформное отображение, несжимаемая жидкость, уравнения Навье — Стокса, условия проскальзывания, вихрь, функция тока.

Библиографическая ссылка: Пивоваров Ю.В. Расчет движения вязкой жидкости, частично заполняющей вращающуюся полость, при больших числах Рей-нольдса // Вычислительные технологии. 2019. Т. 24, № 3. С. 88-105. DOI: 10.25743/ICT.2019.24.3.007.

Введение

Одним из методов очистки полупроводниковых материалов (в частности, германия) является зонная плавка во вращающемся контейнере. Метод состоит в следующем. Цилиндрический контейнер, расположенный под небольшим углом к горизонту, наполнен поликристаллическим полупроводниковым материалом. Контейнер, вращаясь, совершает медленное поступательное движение вдоль своей оси, при этом часть его нагревается до высокой температуры. Материал полупроводника плавится. Образуются фронты плавления и кристаллизации. Расплав не целиком заполняет контейнер: имеется свободная поверхность. Кроме того, его отделяет от стенки контейнера тонкий слой мелкодисперсной смазки.

Опишем схему рассматриваемого процесса. Пусть Oxyz — декартова система координат с началом в центре расплавленной зоны. Ось Oz совпадает с осью контейнера, а ось Оу лежит в вертикальной плоскости, проходящей через ось Oz. Сечение участка контейнера, прилегающего к расплавленной зоне, плоскостью х = 0 показано

* Title translation and abstract in English can be found on page 104.

© ИВТ СО РАН, 2019.

на рис. 1, а, плоскостью ^ = 0 — на рис. 1, б. Здесь Я — радиус контейнера; ^ — угол его наклона к горизонтальной плоскости. Направление поступательного движения совпадает с направлением оси Ох, при этом жидкая зона остается неподвижной, направление вращения контейнера на рис. 1, б —по часовой стрелке.

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

В [1-6] изучались в основном режимы движения, когда жидкость полностью покрывает стенки цилиндра, а внутри имеется газовая полость. Так, в [1, 2] аналитическими методами исследовано поведение слоя жидкости в быстро вращающемся цилиндре, при этом центробежные силы доминируют над капиллярными и гравитационными (действие капиллярных сил не учитывалось). Свободная поверхность жидкости близка к цилиндрической. Экспериментальное изучение такого режима течения проводилось в [3]. В работе [4] при малых скоростях вращения обнаружены трехмерные стационарные течения. Они возникают при увеличении скорости от нуля еще до того, как жидкость полностью покроет боковую поверхность цилиндра. В [5] получены некоторые необходимые условия существования и достаточные условия несуществования плоских течений в двусвязной области в медленно вращающемся цилиндре. В [6] решалась задача о безынерционном движении жидкости во вращающемся цилиндре методом граничных элементов в отсутствие капиллярных сил для случая, когда жидкость полностью покрывает стенки цилиндра.

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

а

б

У 2

3

У

Рис. 1. Сечения контейнера плоскостью х = 0 (а) и плоскостью г = 0 (б): кривые 1-3 соответствуют сечениям фронта кристаллизации, свободной поверхности и фронта плавления

моделировалось действием касательного напряжения, пропорционального разности скоростей стенки и частиц жидкости. Впервые это условие предложено Навье в 1827 г. [9].

К концу 40-х годов XX в. широко использовались условия прилипания. Однако во второй половине XX в. в результате интенсивных экспериментальных исследований неньютоновских жидкостей обнаружено, что для многих из этих жидкостей справедливы соотношения проскальзывания [10]. А в 1982 г. в работе [11] теоретически доказано, что применение граничных условий прилипания даже для ньютоновских жидкостей в случае движущейся линии трехфазного контакта приводит к бесконечному интегралу Дирихле для скорости, если значение угла подхода жидкости к твердой стенке отлично от нуля и значения ж. В связи с этим в работе [12] предложена модель, согласно которой условия прилипания выполняются на всей границе жидкости с твердой стенкой, за исключением некоторых окрестностей линий трехфазного контакта, в которых выполняются условия проскальзывания. В более ранней работе [13] постановка условий проскальзывания в окрестности линий трехфазного контакта при колебательном движении жидкости относительно стенки определялась тем, что жидкость, "смочив в первый период стенки сосуда, при следующих циклах совершает движение по поверхности, покрытой пленкой малой толщины. Это имеет место для частичного или полного смачивания поверхности. Как правило, время стекания пленки существенно больше периода колебаний жидкости". В настоящей задаче пленка в пристеночной области состоит из другого вещества — мелкодисперсной смазки, покрывающей всю поверхность контейнера. Поэтому естественно установить условие проскальзывания на всей границе жидкости с твердой стенкой. Теоретическое обоснование такого подхода дано в [12]. Среди других теоретических работ, посвященных исследованиям задач гидромеханики ньютоновской жидкости при краевых условиях проскальзывания, можно отметить работы Н. Бе1гао аа У^а [14], Н. Еи^а [15], Б. ИоЬ, N. Тапака и А. Таш [16].

Стационарная плоская задача о движении жидкости во вращающемся контейнере решается конечно-разностным методом с использованием разностей против потока при аппроксимации конвективных членов. Для этого аналитически осуществляется конформное отображение полосы на исходную область с использованием биполярной системы координат. Вычисляются первые члены асимптотического разложения вихря и функции тока в окрестности угловых точек области, соответствующих бесконечно удаленным точкам полосы. Затем эта полоса обрезается и на линиях разреза ставятся асимптотические условия для вихря и функции тока, необходимые для замыкания разностной системы уравнений. В [7, 8] в пределе малых скоростей вращения контейнера получены два класса аналитических решений задачи: одно в виде интеграла для специального значения коэффициента проскальзывания, когда задачи для вихря и функции тока расщепляются, и одно в виде ряда для случая, когда область течения является полукругом. Последнее из них используется в настоящей работе для тестирования численного алгоритма.

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

Пусть Г^ = {х,у/у = 0, х Е [— 1,1]} — невозмущенная свободная граница жидкости, Г2 = {х,у/у = — V1 — х2, х Е [-1,1]} — граница жидкости с твердой стенкой контейнера, который вращается по часовой стрелке. Ось х направлена по горизонтали вправо, ось у — вертикально вверх. Область, занимаемую жидкостью, обозначим И. В области И угловые точки х = ±1, у = 0 являются точками скользящего трехфазного контакта.

Отобразим конформно полосу D1 = {£, г]/ — ж < £ < ж, rf < г] < к] на область D посредством преобразования

X

csh£

ch £ — cos rj

с sin Г]

ch — cos

Коэффициент Ламе-отображения равен

Н

ch — cos

В нашем случае, т. е. когда область D является полукругом, имеют место равенства г]° = к/2, с = 1, но для сохранения общности рассуждений будем считать, что с Е (0,1], if = arcsin с. При преобразовании (1) прямая г] = if перейдет в кривую Г2, прямая г/ = к — в отрезок Г1, точка £ = —ж перейдет в левую, точка £ = ж — в правую точки трехфазного контакта. Система координат г/ называется биполярной. Для решения данной задачи она была предложена О.М. Лаврентьевой.

Обозначим через u, v компоненты скорости частиц жидкости в направлениях г/. Введем вихрь ш и функцию тока Ф следующим образом:

_ id (Ну ) д (Ни) ш = Н2 V ~д£

_ 1 дФ Н d

1 дФ

" = — НЖ

Задача для определения этих функций в биполярной системе координат [7, 8]:

ДФ = —Н 2и,

„ (дФдш дФдш \ _ . ^

Лш = ЧаЦаЦ — ЖаШ,)' К^ Е1,1-

ф(е.ч0) = о, ФК,к) = о, ш('''0) = —A1 + (2 — Al)Н«^ дфФ

ш(е ,к) = о,

2 2 где А = — +

д^2 дг]2

Задача для определения формы возмущенной свободной границы: Пх) — Б1о/"(Ж) = Ф(Ж) + /с, X Е (—с, с),

(4)

(5)

(6) (7)

№)U=-e = — ^, f'(x)lx=c = М ,

f(x)dx = 0,

:iq)

где

Ф(х(0)

Fr

У

1 дФ Н

Fr Re

dr] ^ Нд£[ тт

(_1 дФ ^

Г1=Ж

с

2

Здесь уравнение (3) — следствие (2); (4) — уравнение импульса; (5) — условия непротекания на твердой стенке и свободной границе; (6) — условие проскальзывания Навье (касательные напряжения пропорциональны разности скоростей частиц жидкости и стенки контейнера с коэффициентом пропорциональности, равным —к); (7) — условие отсутствия касательных напряжений на свободной границе; (8) — условие капиллярного равновесия; (9) — условие подхода свободной границы жидкости к твердым стенкам под углом п/2; (10) — закон сохранения массы жидкости; /(ж) — функция, описывающая форму возмущенной свободной границы жидкости, /0 — константа, подлежащая определению.

Безразмерными критериями подобия являются

Ще = ™, Рг = Ш, Во = А1 = (11)

V д0К а ри

где и — скорость движения стенок контейнера, К = 1 см — радиус контейнера, V = 1.35е — 3 см2/с — кинематическая вязкость жидкости (германий), д0 = 9.81 см/с2 — ускорение силы тяжести, р = 5.571 г/см3 — плотность жидкости (германий), а = 600 дин/см — поверхностное натяжение жидкости (германий), к — коэффициент проскальзывания, И,е — число Рейнольдса, Рг — число Фруда, Во — число Бонда, А1 — параметр проскальзывания.

В условиях эксперимента скорость движения стенок контейнера составляла примерно 60 см/с (10 оборотов в секунду), что приводит к числу Рейнольдса порядка 45 000 и числу Фруда Рг = 3.6. Но характерная скорость движения самой жидкости была примерно в 10 раз меньше. Обозначим через vmax отношение максимума скорости движения расплава и скорости движения стенок контейнера. Тогда реальное число Рейнольдса будет И,ед = утахИ,е = 4500, а реальное число Фруда — Ргд = ^ахРг = 0.036. Число Бонда Во = 9.1, параметр проскальзывания А1 положим равным 0.2.

2. Построение последовательности ортогональных разностных сеток

Пусть П = ,г\/ — £0 <£< г/0 < г/ < п} — прямоугольная область, полученная из области посредством отбрасывания бесконечных ее участков справа и слева от разрезов £ = ±£0. При расчетах будем полагать £0 = 9.

При построении сетки в П необходимо сгущать ее в окрестности линии £ = 0, г/ Е [ г]0, п], чтобы получить более или менее равномерную сетку в области И2, являющейся образом прямоугольника П при преобразовании (1). Кроме того, ее необходимо сгущать в окрестности отрезков г] = гр и г] = п, ^ Е [—£0,£0], так как в окрестности линий Гх и Г2 образуются пограничные слои, связанные с большим числом Рейнольдса задачи.

В работе [17] описан алгоритм построения одномерной сетки, сгущающейся на границах отрезка [0,1]. Для его применения необходимо задать шаги сетки на левом и правом концах отрезка и число разбиений. Построим согласно этому алгоритму сетку ап, п = 0, N/2, сгущающуюся около правой границы отрезка [0,1] так, что

«1 = 2hal/N, ам/2 — ам/2-1 = 2^2 /N.

Далее полагаем

ап :=ап/2, п = 0,^2, ап := 1 — ам_га, п = N/2 + 1^.

В результате получится сетка, сгущающаяся в середине отрезка [0,1]. После этого вычислим узлы

:= - 1), п = о;ж

Затем построим сетку , п = 0,М, сгущающуюся на границах отрезка [0,1], такую, что

А = /м, рм - рм_ = ^ /м,

и вычислим узлы

цт := г]° + (п - 'ц°)рт, т = 0,М.

Числовые значения управляющих параметров при построении описанных выше сеток приведены ниже:

щ Щ Л? Л?

1 0.25 0.325 0.5

Далее строится последовательность сеток следующим образом. Сначала определяются узлы сеток при N х М = 32 х 8. Затем число узлов то по направлению 'ц, то по направлению £ удваивается и строятся узлы новых сеток до тех пор, пока не будет построена сетка размерности 2048 х 512.

На рис. 2 в качестве примера приведена сетка размерности 128 х 32 в области , являющаяся образом сетки в прямоугольнике П при замене переменных, определяемой формулами (1).

3. Построение асимптотики вихря и функции тока вблизи точек контакта

При выводе асимптотических формул целесообразно сначала рассмотреть случай, когда гр = ^/2 + е, а затем перейти к пределу при е ^ 0. Требуется определить поведение функций Ф, ш при |£| ^ то. Старший член асимптотики для вихря удовлетворяет задаче

Аш° = 0, (С/п) е £>ь ш°(С,'П°) = -Л!, ш°(£,*) = 0.

Ее решением является функция

ж — 'П

ш°(£,ч) = ш°(ц) = -Л-°. (12)

ж -

Подставляя ее в правую часть (9) и учитывая, что Н2(£, г]) ~ 4с?е 2|^, получим асимптотическую задачу для Ф:

= 4А11—Л .-21?1

ДФ „ = 4 с2 А1

Ж — Г]с

(С, V) ЕБи

Ф0(£, г]с) = 0, Ф0(£,^) = 0.

Возьмем частное решение неоднородной задачи Фя и набор решений однородной задачи Ф^, причем каждое решение терпит разрыв производной при £ = 0. Составим их линейную комбинацию так, чтобы вся сумма не имела этого разрыва:

Фя = С — е-2|С|,

Н \ж — Г]с вт(2 г]с)) '

Ф° = 81п( к 3 С(Г] — Г]с)) е—к^,

те

Ф« = Ф^ — ^ акФ°к,

к=С

где 30

Ж

Ж — Г]

. Потребуем, чтобы

д Ф,

д

?=С

Получим

с2А1

(

те

кж

Ж — Г] 81п(2г]) \ ж — Г]с — 81п(2г7С^ = ^ ак 2(тг — Г]с)

яп( к 3с (Г1 — 'Пс)).

Разложим в ряд Фурье функцию

ж — г]с 8т(2?7С) Величины Ьк оказываются равными

^с — —72^ = ЕЬк 81п(^3 ^ — г1)).

к= 1

8

жк(4 — (3С )2к2)

13)

14)

тогда

Итак, получено

2 кл 2Ьк ак = с А1

16 с2А1

Ф«(£, V)

3Ск ж3С к2(4 — (3С )2к2)'

К — У _ е-2|?|_

ж — г]с вт^?/0) у

Е

81п( к3°(Г] — Г]°)) | С2А1 = Фв + 0(е—2^0|?|),

16

ж3С ^ к2(4—(3С)2к2)

15)

0

к

, V) - rqo sin(2rq0) J nß0(4 - (ß0)2) e j c AL (1б)

Пусть теперь

""=2+- ß0=ên ■ (17)

тогда sin(2^0) ^ 0, 4 — (ß0)2 ^ 0. Переходя в (15) к пределу при е ^ 0, запишем

lim Ф a = Фа + 0(е-4|ç|),

c2Ale-2IÇI

Фа =-¡т-{(5 + 4|£|) sin(2rç) + 4(тт — rç)(1 + 008(2/7))} . (18)

2n

Найдем асимптотику скорости на границе г] = г]0 при е > 0:

e|ç| <9Фa

Va(0

2с drq

va(0 + 0(e

«■«> =**+• (19)

Асимптотическая задача для следующего члена разложения функции ш имеет вид

Д£0 = 0, (£, V) ша(£, г,°) = (2 - Л1)0в(0, £«(£,*) = 0.

Решив ее, получим

¿а = ¿а + 0(е-ß"|£| ),

= — й—А»^ ((_!_ + 2ctg(24»)) e-и +

2 — Г]0 ) 8in(n — Г]0)

+ 16 sin((ß0 — 1)(n — y)) e-(ß0-1)^ (20)

+ n sin((ß0 — 1)(тг — r]0)) (4 — (ß0)2)' . (Uj

1

Переходя в (19), (20) к пределу при e ^ 0, с учетом (17) найдем, что V ~ (5 + 4|С|)е-|ç| .. .. _

lim Va(0 = Va(0, Va(£) =--"-cAl, lim Ша = ¿a,

s—^0 2n s—0

Ша = — (2 — Al)cAl {(5 + 4|£|) sin(n — ri) — 4(тг — Г]) cos(n — Г])} e-|ç|. (21)

2n

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

. Ф /<9Ф адш 0 д Ф а дщ>\ )сЛ (22)

= ^Г— Wj , (е,v) eDl, (22)

n=n0

Шя(£, Vе) = 0, Шп(С,п) = 0.

(23)

Здесь Фа определяется формулой (15), ш- — формулой (12). Используя разложение (13) с учетом (14), представим правую часть (22) в виде

д Ф 0 дш- 16 c2A12sign(e)Re ^ е-^-е-^ . - -

д д

п(тг — Г]-) ^ к((3-)2к2 — 4)

Решение задачи (22), (23) имеет вид

те ж

— те п0

ип = / "¡Са,п,1 /) (—ЯедФЬ

где

? 1 1п еЬ(3-(| — 0) — сс8(3-(у + у — 2г]-))

, п еЬ(3-(| — 0) — сов(3-('Ц — V))

те ,—к/З01?—

£

к= 1

кж

яп( к3-(г] — Г]-)) яп( к3-(г1 — Г]-))

функция Грина для уравнения Лапласа в И1. Используя формулы

в1п2( к 3-(г] — г]-))((г] =

Ж — Г]-2 ,

яп(к3-(г] — Г]-)) яп(п3-(Г] — г]-))((г] = 0, к = п

—к130 —

=

4е—fc^0|g| — 2 к3-е—2 4 — к2(3-)2

—те те

е—к/З01?—?|е—к?0

получим

Шп =

4 с £

к=1

4 е— 2 к 3-е—2|€| к2(3-)2 — 4

+ (ж- + и)

э1п( кЗ- (Г] — Г]-)) к2((3-)2к2 — 4)

2 c2A12sign(£)Re

С = -

ж2

(24)

ж

ж

Отсюда следует, что

шн = ШК + 0(£ е—2/°^)Яе,

П 8т((р°(г! - г]0))

Шк = 4С

4 е—/°|£|_2 В°е-2|«1 /1 \ ,

(/°)2 - 4 + + К1)е

( / °)2 - 4

о д°

2/° в1п( к/°('п - т)°))

к((/°)2к2 - 4)2

е-21?1

Переходя к пределу при е ^ 0, получим

Иш шя = шК + 0(С е-4|?|)Яе,

°

шя = С {8 (4|С!2 + 2|£| + 1) йп(ж - 2Ч) - £ 8'П(к(2-" } ^^

(25)

Заметим, что ряд в (25) легко посчитать численно, так как его члены убывают со скоростью 0(к-5).

Таким образом, построение асимптотического решения для вихря и функции тока закончено.

При расчетах на линиях разреза полосы £ = ±£° будем ставить условия первого рода

Ф(±£°, ??) = Фа(£°, V), (26)

ш(±£°, V) = ш°('П) + ша(£°, п) + шК(±£°, п), (27)

где функция Фа определена формулой (18), а функции ш°, ша, — соответственно формулами (12), (21), (25) с учетом (24).

4. Тестовые расчеты

При г]° = ж/2, И,е = 0 рассматриваемая задача имеет точное решение [7, 8]

ФА = А1 г, ф)(1 - г2), ша = 2А1(2 - А1)^(г, ф) + г, ф),

N 2 ^ г2™-1 яп((2т - 1)ф) , Л

= 2 £ (2т - 1)(4'„ + А1 -ф2) ■ <28>

т=1 4 ' 4 '

2А1 ( 2 гвтф \ 2А1(ж - г])

г-г / \ (2г

^2(^ ф) = — аг^ I 2 I

ж \ 1 - г2 / ж

г = л/х2 + у2, ф = - агссов(х/г), вЬ £ вт г]

х = Т7-, у

сЬ £ - сов г] сЬ £ - сов г]

Число членов ряда Мтах в (28) определялось по формуле

М Г -[3/^(г)] при г < 0.9998, тах = \ -[3/ 1ё(о.9998)] при г > 0.9998.

Квадратные скобки означают взятие целой части.

Основные расчеты произведены при асимптотических условиях (26), (27). При использовании граничных асимптотических условий третьего рода [7, 8] в задаче для функции тока сходимость итераций была значительно хуже, а точность решения — та же.

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

еФ = // 1Ф - |Н2^, £ш = Ц -^|Н2^.

п п

В табл. 1 приведены значения величин еф, еш в зависимости от чисел разбиений N по оси £ и М по оси г], а также порядки сходимости рф, рш, определяемые по формулам

Рф = ^2(£ф2/£Ф1), Рш = ^2(^2 /еШ1), (29)

где £ф2, £ф1 — значения £ф при N = N1, М = М1 и N = 2Nl, М = 2М1 соответственно. Значения £ш2, £ш1 определялись аналогично. Так как порядки сходимости при изменении чисел разбиений сильно меняются, имеет смысл определить средние по всей совокупности чисел разбиений порядки сходимости

5 5

РФ = ^ РФг/5, рш = ^ Р^г/5, г=1 г=1

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

Таблица 1. Порядки сходимости функции тока и вихря к точному решению при Ив = 0

г N х М £ф £ш PФi Ршг

32 х 16 3.82е - 4 9.54е - 3 — —

1 64 х 32 2.58е - 4 5.56е - 3 0.556 0.779

2 128 64 1.45е - 4 2.99е - 3 0.831 0.895

3 256 х 128 2.55е - 5 5.52е - 4 2.507 2.437

4 512 х 256 1.65е - 5 3.07е - 4 0.628 0.846

5 1024 512 1.17е - 5 1.78е - 4 0.496 0.786

Среднее — — 1.003 1.149

Таблица 2. Порядки сходимости "в себе" функции тока и вихря при Ив = 0

г N х 2М £ф £ш PФi Ршг

32 16 4.28е - 4 4.30е - 3 — —

1 64 32 2.62е - 4 4.34е - 3 0.708 -0.01

2 128 х 64 1.46е - 4 2.69е - 3 0.846 0.691

3 256 128 1.30е - 4 2.48е - 3 0.173 0.114

4 512 256 1.25е - 5 2.49е - 4 3.377 3.316

5 1024 х 512 5.02е - 6 1.27е - 4 1.312 0.977

Среднее — — 1.283 1.018

Таблица 3. Порядки сходимости "в себе" функции тока и вихря в зависимости от числа Рейнольдса

Ив £ф £ш Рф рш

10 1.39е - 5 2.76е - 4 — —

5.47е - 6 1.39е - 4 1.345 0.993

1000 6.02е - 5 1.18е - 3

3.07е - 5 6.07е - 4 0.974 0.956

6000 2.01е - 4 3.86е - 3

1.10е - 4 2.16е - 3 0.871 0.839

20 000 3.23е - 4 6.43е - 3

1.85е - 4 3.73е - 3 0.801 0.783

45 000 4.02е - 4 8.15е - 3

2.39е - 4 4.84е - 3 0.752 0.753

Таблица 4. Порядки сходимости "в себе" функций Ф, / в зависимости от числа Рей-нольдса

Ив £ф Рф Pf

10 5.87е - 4 2.91е - 3 — —

2.92е - 4 1.44е - 3 1.006 1.013

1000 7.49е - 3 9.36е - 3

3.81е - 3 4.74е - 3 0.973 0.982

6000 2.56е - 2 2.65е - 2

1.41е - 2 1.46е - 2 0.862 0.866

20 000 4.56е - 2 4.69е - 2

2.56е - 2 2.63е - 2 0.832 0.833

45 000 5.96е - 2 6.10е - 2

3.36е - 2 3.45е - 2 0.828 0.825

В табл. 2 приведены порядки сходимости "в себе", определяемые также по формулам (29), в которых величины е^2, £ф1, еш2, £ш1 имеют другой смысл. А именно, есть интегральная норма разности между двумя приближенными решениями, построенными при N = М1, М = М1 и N = 2М1, М = 2М1 (сначала массивы функции тока размерностей N1 х М1 и 2N1 х 2М1 приводятся к размерности N1 х 2М1, что позволяет использовать только одномерную интерполяцию, а затем сравниваются), а еф1 — аналогичная величина для N = 2N1, М = 2М1 и N = 4^, М = 4М1. Величины еш2, еш1 определяются аналогично.

Как видно из табл. 1,2, порядки сходимости к точному решению и порядки сходимости "в себе" при И,е = 0 в среднем очень близки к единице.

В табл. 3 приведены интегральные нормы разности еФ2, еФ1, еш2, еш1 между двумя приближенными решениями, построенные, как это описывалось выше, при N1 = 512, М1 = 128, и порядки сходимости "в себе" рф, рш в зависимости от числа Рейнольдса. Из табл. 3 видно, что величины рф, рш монотонно убывают с увеличением числа Рейнольдса, при этом их отклонение от единицы в меньшую сторону составляет не более одной четверти.

В табл. 4 приведены относительные интегральные нормы разности еф, еf между

двумя приближенными решениями, где еф определялась по формуле:

1 1

£Ф = J |Ф1(х) - Ф2(х)^х/ J |Ф1(х)|^х, -1 -1

а еf — по аналогичной формуле с заменой Ф на /, Ф ^ — правые части уравнения (8), /г(х) — решения задачи (8)-(10), а также порядки сходимости "в себе" для этих функций при N1 х 2М1 = 512 х 256. Из табл. 4 видно, что все порядки сходимости не сильно отличаются от единицы.

5. Результаты расчета гидродинамического течения

Задача решалась методом конечных разностей. При расчете функции вихря в левую часть уравнения (4) вводился член дш/дЬ и решение стационарной задачи находилось методом установления по времени. Для решения эволюционной задачи использовался

метод приближенной факторизации. При аппроксимации конвективных членов в уравнении импульса использовались разности против потока, обеспечивающие монотонность схемы при любом шаге по времени за счет диагонального преобладания в матрицах системы [18]. Применялся алгоритм, позволяющий точным образом выполнить условие (6) на каждом шаге по времени [19]. При решении задачи для функции тока использовался метод В.Г. Зверева [19, 20]. Задача определения свободной границы решалась методом функции Грина, а входящие в решение интегралы от функции Ф(ж) рассчитывались численно [7, 8].

Расчет проводился на последовательности сеток. Сначала функции Ф, ш рассчитывались на сетке 32 х 8. При этом в качестве начального условия для функции тока использовалась функция Фа, для вихря — функция шо + ша. На первых двухстах итерациях по времени число Рейнольдса линейно увеличивалось от нуля до заданного значения, а затем оставалось постоянным. Затем решалась задача на сетке 32 х 16 с начальным приближением, полученным путем интерполяции с ранее построенного решения на сетке 32 х 8. Процесс удвоения числа узлов сетки то по одному, то по другому направлению с последующим нахождением функций Ф, ш продолжался до тех пор, пока не было построено решение на сетке размерности 2048 х 512.

Шаг по времени для всех вариантов был равен 1.0 е - 4. Число итераций по времени для максимального числа разбиений и максимального числа Рейнольдса составляло 2000. Время счета — одни сутки. Для всех вариантов максимум модуля производной по времени от функции вихря в конце расчетов составлял не более 1. 0 - 4. При решении задачи для функции тока после каждой итерации вычислялась норма невязки. Затем, если ее значение было меньше 2.0 е - 4 или число итераций Ь для функции тока становилось равным десяти, итерации прекращались. Во всех вариантах к концу расчетов величина Ь была равна единице.

В табл. 5 приведены максимумы модуля безразмерной скорости утах, реальное число Рейнольдса И,ед = ^тахКе, формальное число Фруда Ег (см. формулы (11)), реальное число Фрула Егд = ь^ахЕг, максимум модуля функции тока Фтах, максимум модуля функции вихря штах и максимум модуля отклонения функции свободной поверхности от нуля /тах в зависимости от формального числа Рейнольдса И,е, определяемого по первой формуле (11).

На рис. 3 приведены линии тока при числах Рейнольдса от 10 до 6000. При больших числах Рейнольдса линии тока почти не меняются. На рис. 4 показаны изолинии функции вихря во всем диапазоне чисел Рейнольдса. Видно, что при достаточно больших числах Рейнольдса в окрестности границы области формируется вихревой пограничный слой, затем идет переходная область и далее — область, где значения функции вихря близки к константе.

Таблица 5. Характеристики течения при различных числах Рейнольдса

Ив ^тах Ивд Рг Ргд Фтах ^тах /тах

0 0.1053 0 0 0 0.0215 0.3895 0

10 0.1052 1.053 1.86е - 7 2.06е - 9 0.0215 0.3894 1.14е - 8

1000 0.0836 83.62 1.86е - 3 1.30е - 5 0.0186 0.3505 2.34е - 6

6000 0.0728 436.7 6.69е - 2 3.54е - 4 0.0178 0.3282 6.00е - 5

20 000 0.0735 1470 0.7439 4.02е - 3 0.0175 0.3218 6.84е - 4

45 000 0.0738 3321 3.7659 2.05е - 2 0.0174 0.3193 3.41е - 3

а б в

а б в

а б в

Ив = 45 000 (д)

На рис. 5 показаны функции /(х), описывающие форму свободной границы при числах Рейнольдса от 10 до 6000. При больших числах Рейнольдса форма свободной границы почти не меняется, меняется только ее масштаб по вертикальной оси. Значения максимума модуля /тах функции свободной границы для всех вариантов чисел Рейнольдса приведены в табл. 5.

б.Ое - 5

Рис. 5. Форма свободной границы при Ив = 10 (а), Ив = 1000 (б), Ив = 6000 (в)

Следует отметить, что для каждой из построенных сеток решение стационарной задачи ищется методом установления по времени, так что начальное приближение для функций Ф, ш есть Ф8(х, у) + ДФ(х, у), ш3(х, у) + Дш(х, у), где Ф5(х, у), ш3(х, у) — искомое стационарное решение, а ДФ(х, у), Дш(х, у) — его возмущение. Установлено, что изолинии функций ДФ(х, у), Дш(х, у) повторяют изолинии функций Ф5(х, у), ш3(х, у), но имеют значительно меньшую амплитуду. Так, например, при N = 2048, М = 512, И,е = 45 000 возмущения функций Ф, ш имеют амплитуды 1.63е - 4 и 3.18е - 3 соответственно при амплитудах функций Ф5, ш3, равных 1.74е - 2 и 3.19е - 1 (табл. 5). Это низкочастотные возмущения. Кроме того, есть еще высокочастотные возмущения, связанные с ошибками интерполяции при увеличении числа разбиений области по одному из направлений. Таким образом, можно говорить об устойчивости стационарного решения к описанным выше возмущениям. Вопрос об устойчивости решения к другим типам возмущений остается открытым.

Заключение

Решена задача о плоскопараллельном стационарном движении вязкой несжимаемой жидкости, частично заполняющей цилиндрическую вращающуюся полость, для случая, когда область течения с невозмущенной свободной границей является полукругом. При решении использованы метод конформного отображения при построении ортогональной разностной сетки, метод установления по времени при решении стационарной задачи, метод приближенной факторизации при решении эволюционного уравнения для вихря, метод В. Г. Зверева при решении задачи для функции тока, метод, позволяющий точно выполнить граничные условия, связывающие вихрь на границе и приграничные значения функции тока на каждом шаге по времени, и конечно-разностный метод. Расчет проводился на последовательности ортогональных разностных сеток, начиная с сетки размерности 32 х 8 и кончая сеткой 2048 х 512. Построены линии тока, изолинии функции вихря течения, а также форма возмущенной свободной границы. Получено, что при больших числах Рейнольдса область течения делится на три части: вихревой пограничный слой в окрестности границы, переходную зону и зону, в которой функция вихря имеет постоянное значение. Появление переходной зоны обусловлено, по-видимому, наличием угловых точек границы области, а также тем, что в этих точках функция вихря терпит скачок.

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

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

[1] Ждан С.А. Задача о движении вязкой жидкости во вращающемся круге в поле тяжести // Вестн. Моск. ун-та. Сер. 1. Математика, механика. 1987. № 1. С. 86-89.

Zhdan, S.A. The problem of movement of viscous fluid in a rotating circle in the field of gravity // Vestn. Mosk. Un-ta. Ser. 1. Matematika, Mekhanika. 1987. No. 1. P. 86-89. (In Russ.)

[2] Greenspan, H.P. On a rotational flow disturbed by gravity // J. of Fluid Mech. 1976. Vol. 74. Pt 2. P. 335-352.

[3] Gans, R.F., Yalisove, S.M. Observations and measurments of flow in partially-filled horizontally rotating cylinder // Trans. ASME. Ser. 1. J. of Fluids Eng. 1982. Vol. 104, No. 3. P. 363-366.

[4] Balmer, R.T., Wang, T.G. An experimental study of internal hidrocyts // Trans. ASME. Ser. 1. J. of Fluids Eng. 1976. Vol. 98, No. 4. P. 688-694.

[5] Бадратинова Л.Г. О движении жидкого слоя по внутренней поверхности горизонтального вращающегося цилиндра // Динамика сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. Ин-т гидродинамики. 1993. Вып. 106. С. 179-184.

Badratinova, L.G. On motion of fluid layer over an internal surface of horizontal rotating cylinder // Dinamika Sploshnoy Sredy. 1993. Iss. 106. P. 179-184. (In Russ.)

[6] Шрагер Г.Р., Козлобродов А.Н., Якутенок В.А. Моделирование гидродинамических процессов в технологии переработки полимерных материалов. Томск: Изд-во Томск. ун-та, 1999. 230 с.

Shrager, G.R., Kozloborodov, A.N., Yakutenok, V.A. Simulation of hydrodynamic processes in technology of refinement of polymeric materials. Tomsk: Izd-vo Tomsk. Un-ta, 1999. 230 p. (In Russ.)

[7] Проблемы вычислительной математики / А.Ф. Воеводин, В.В. Остапенко, Ю.В. Пивоваров, С.М. Шугрин. Новосибирск: Изд-во СО РАН, 1995. 154 с.

Problems of computational mathematics / A.F. Voevodin, V.V. Ostapenko, Yu.V. Pivovarov, S.M. Shugrin. Novosibirsk: Izd-vo SO RAN, 1995. 154 p. (In Russ.)

[8] Пивоваров Ю.В. Моделирование конвекции расплава полупроводникового материала при зонной плавке: Дис.... канд. физ.-мат. наук. Новосибирск, ИГиЛ СО РАН. 2006. 135 с. Pivovarov, Yu.V. Convection modelling of a semiconductor material in melting zone: Dis. ... kand. fiz.-mat. nauk. Novosibirsk, IGiL SO RAN. 2006. 135 p. (In Russ.)

[9] Navier, C.-L.-M.-H. Memoire sur les lois de l'equilibre et du mouvement des corps solides elastiques (1821) // Mem. Acad. Sci. Inst. France. 1827. Vol. 7(2). P. 375-393.

[10] Астарита Дж., Марручи Дж. Основы гидромеханики неньютоновских жидкостей. М.: Мир, 1978. 309 с.

Astarita, G., Marruchi, G. Principles of non-Newtonian fluid mechanics. McGraw - Hill, 1974. 296 p.

[11] Пухначев В.В., Солонников В.А. К вопросу о динамическом краевом угле // Прикл. матем. и механика. 1982. Т. 46, № 6. С. 961-971.

Pukhnachev, V.V., Solonnikov, V.A. On the problem of dynamic contact angle // J. of Appl. Math. and Mechanics. 1982. Vol. 46, No. 6. P. 771-779.

[12] Байокки К., Пухначев В.В. Задачи с односторонними ограничениями для уравнений Навье —Стокса и проблема динамического краевого угла // ПМТФ. 1990. Т. 31, № 2. С. 27-40.

Baiocchi, C., Pukhnachev, V.V. Problems with one-sided constraints for Navier — Stokes equations and the dynamic contact-angle problem // J. of Appl. Mechanics and Techn. Physics. 1990. Vol. 31, No. 2. P. 185-197.

[13] Богоряд И.Б. Динамика вязкой жидкости со свободной поверхностью. Томск: Изд-во Томского ун-та, 1980. 102 с.

Bogoryad, I.B. Dynamics of viscous fluid with a free boundary. Tomsk: Izd-vo Tomsk. Un-ta, 1980. 102 p. (In Russ.)

[14] Beirao da Veiga, Н. Regularity for Stokes and generalized Stokes systems under nonhomogeneous slip-type boundary conditions // Advances in Differential Equations. 2004. Vol. 9, No. 9-10. P. 1079-1114.

[15] Fugita, H. Non-stationary Stokes flows under leak boundary conditions of friction type // J. of Comput. Mathematics. 2001. Vol. 19. P. 1-8.

[16] Itoh, S., Tanaka, N., Tani, A. The initial value problem for the Navier-Stokes equations with general slip boundary conditions in Holder spaces // J. of Math. Fluid Mechanics. 2003. Vol. 5. P. 275-301.

[17] Пивоваров Ю.В. Расчет процесса сближения двух сферических капель, находящихся в среде Бингама // ПМТФ. 2014. Т. 55, № 5. С. 29-44.

Pivovarov, Yu.V. Calculating the approach of two spherical droplets located in a Bingham fluid // J. of Appl. Mechanics and Techn. Physics. 2014. Vol. 55, No. 5. P. 750-763.

[18] Пивоваров Ю.В. Условия монотонности факторизованной разностной схемы для эволюционного уравнения с двумя пространственными переменными // Вычисл. технологии. 2001. Т. 6, № 4. С. 81-91.

Pivovarov, Yu.V. Monotonicity conditions of factored difference scheme for an evolutionary equation with two space variables // Comput. Technologies. 2001. Vol. 6, No. 4. P. 81-91. (In Russ.)

[19] Пивоваров Ю.В. Расчет движения жидкости с переменной вязкостью в области с криволинейной границей // Вычисл. технологии. 2005. Т. 10, № 3. С. 87-107. Pivovarov, Yu.V. Calculation of a flow of liquid with variable viscosity in a domain with curvilinear boundary // Comput. Technologies. 2005. Vol. 10, No. 3. P. 87-107. (In Russ.)

[20] Зверев В.Г. Об одном итерационном алгоритме решения разностных эллиптических уравнений // Вычисл. технологии. 1999. Т. 4, № 1. С. 55-65.

Zverev, V.G. An iterative algorithm for the solution of difference elliptical equations // Comput. Technologies. 1999. Vol. 4, No. 1. P. 55-65. (In Russ.)

Поступила в 'редакцию 3 октября 2018 г., с доработки — 8 ноября 2018 г.

Computation of motion of viscous fluid partially filling a rotating cavity at large Reynolds numbers

PIVOVAROV, YURIY V.

Lavrentyev Institute of Hydrodynamics SB RAS, Novosibirsk, 630090, Russia Corresponding author: Pivovarov, Yuriy V., e-mail: pivov@hydro.nsc.ru

The goal of the research is to simulate a hydrodynamic flow of semiconductor material melt during its purification from impurity by using the method of zone melting in a horizontal cylindrical rotating container. The container is not fully filled with melt, there is a free surface. The flow is considered to be plane-parallel and stationary. It is supposed that the flow domain with unperturbed free boundary is a semicircle. The melt is separated from a container wall by a thin layer of finely dispersed lubricant. Therefore, when stating a problem, the Navier slip condition is set on the boundary with a wall.

The problem formulation in variables of vortex-stream function is performed. To solve it, the following methods are used, namely, relaxation method for time, the method of an approximate factorization for solution of an evolutionary equation for vortex,

© ICT SB RAS, 2019

V.G. Zverev's method is used to solve the problems for stream function, this method allows precise satisfying a boundary condition which connects a vortex on boundary and near-boundary values of stream function at each time step using finite-difference method and method of computation on the sequence of meshes, beginning with the mesh of dimension 32 x 8 and ending with the mesh of dimension 2048 x 512.

The research results in constructed patterns of streamlines and vortex isolines along with the form of perturbed free boundary at various values of Reynolds number. The range of Reynolds number is from zero to three thousand that corresponds to the experimental data.

Conclusions.

1. The problem on motion of viscous fluid filling half the cylindrical rotating cavity is solved.

2. A great number of mesh points made it possible to minimize the influence of scheme viscosity and to obtain trustworthy results at large Reynolds numbers.

3. It is shown that the flow domain is divided into a vortex boundary layer, transition zone and the zone of constant vortex.

Keywords: conformal mapping, incompressible fluid, the Navier —Stokes equations, slip condition, vortex, stream function.

Cite: Pivovarov, Yu.V. Computation of motion of viscous fluid partially filling a rotating cavity at large Reynolds numbers // Computational Technologies. 2019. Vol. 24, No. 3. P. 88-105. (In Russ.) DOI: 10.25743/ICT.2019.24.3.007.

Received October 3, 2018 Received in revised form November 8, 2018