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

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

Том 17, № 4, 2012

Расчёт движения осесимметричного твёрдого тела

в жидкости Бингама*

Ю. В. ПИВОВАРОВ Институт гидродинамики СО РАН, Новосибирск, Россия e-mail: [email protected]

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

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

Введение

Рассматривается нестационарное движение осесимметричного твёрдого тела в жидкости Бингама (матрице) вдоль оси симметрии под действием силы тяжести. Возникающее течение в матрице также предполагается осесимметричным. Хорошо известен тот факт, что вне некоторой поверхности вокруг рассматриваемого тела поле скоростей в матрице будет тождественно равно нулю. Поэтому можно рассматривать движение жидкости внутри такой поверхности, например, шара достаточно большого радиуса, эллипсоида и т. д. В цилиндрических координатах r, z, p, связанных с телом, где r — полярный радиус, p — полярный угол, от которого искомые функции не зависят, ось z направлена вдоль оси симметрии, можно решать задачу в полуплоскости r > 0, p = const. Тогда двумерная область в координатах r, z между охватывающей тело поверхностью, вне которой жидкость покоится, и границей тела будет криволинейным четырёхугольником. Известно, что криволинейный четырёхугольник можно единственным образом конформно отобразить на прямоугольник с заранее неизвестным отношением сторон. В некоторых частных случаях это удаётся сделать аналитически, а в общем случае есть численный алгоритм, основанный на минимизации некоторого функционала [1, 2]. Итак, пусть данное отображение построено. Далее модель жидкости Бингама регуляризуется так же, как в работе [3], где доказано, что при параметре регуляризации 8, стремящемся к нулю, решение регуляризованной задачи сходится к решению исходной задачи для жидкости Бингама. В сущности регуляризованная модель сводится к описанию движения неньютоновской жидкости, у которой вязкость определённым образом зависит от второго инварианта тензора скоростей деформаций. Решение данной задачи

* Работа выполнена при поддержке программы ОЭММПУ РАН (проект № 2.13.2) и РФФИ (грант № 12-01-00149-а).

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

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

1. Размерные параметры задачи

Пусть щ — кинематическая вязкость материала матрицы, м2/с, ко — ее предел текучести, Па, тоо — радиус шара, имеющего такой же объем, что и рассматриваемое тело, м, д = 9.81 м/с2 — ускорение свободного падения, р1 — плотность тела, кг/м3, р0 — плотность материала матрицы. Характерную скорость определим по формуле Стокса [5]

лг _ 2 дгт(р1- ро)

У0 = А -,

9 щро

а характерное время — по формуле

*о = ТтМ.

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

VI = коТоо/(роУо),

а в качестве характерного напряжения — величину ко. Тогда характерная сила определяется как интеграл по поверхности сферы от касательного напряжения, равного ко:

Со = п2Г)око.

2. Безразмерные критерии подобия

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

т, Voroo ЭпМо (pi - po)tog . Voto vo Re =-, Ai = --—, a2 =-T-,-, Ao = -t- , В =-, a =—.

vi 4piroo Vo piVo Ai roo vi

Здесь Re — число Рейнольдса, Ai — отношение характерной силы к произведению массы тела на его ускорение, A2 — отношение разности силы тяжести и силы плавучести к произведению массы тела на его ускорение, Ao — отношение разности силы тяжести и силы плавучести к характерной силе, определяемой пределом текучести материала матрицы, В — кинематический параметр, a — обратное число Бингама [6].

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

Пусть г, г, р — система координат, о которой говорилось во введении. Форма тела описывается параметрически:

г = г^х), г = ¿1 (х), х Е [0,Х]. (1)

Внешняя граница области течения также определяется параметрически:

г = г2(х), г = (х), х Е [0,X]. (2)

Заметим, что поскольку система координат г, г, р движется вместе с телом, то в формулах (1), (2) нет зависимости от времени. Функции

г = гз(х,у), г = ¿з(х,у), х Е [0,Х], у Е [0,У]

конформно отображают прямоугольник 0 < х < X, 0 < у < У на область между границами (1) и (2) так, что

гз(х, 0) = г1(х), гз(х, 0) = ^(х), гз(х,У) = г2(х), гз(х,У) = ^(х), х Е [0,Х],

гз(0,у) = 0, гз(X,у) = 0, у Е [0,У].

Функция безразмерной вязкости регуляризованной модели жидкости Бингама имеет вид [3]

( а + 1/(^ВД) при |£| > 35/2, V(|Я|)= < /(|£|) при 5/2 < |£| < 35/2, (3)

[ а + 1/(^25) при |£| <5/2.

Здесь |Д|2 — второй инвариант тензора скоростей деформаций, 5 — параметр регуляризации (малое число),

/(|Я|) = а + 1/(^25) - 17^2(|Я| - 5/2)з/(2754)+

+35^2(|£| - 5/2)4/(5455) - 5^2(|Я| - 5/2)5/(2756) —

полином пятой степени, осуществляющий склейку значений функции V(|^|), ее первых и вторых производных в точках |Д| = 5/2 и |Д| = 35/2. Граница жидкой зоны определяется соотношением

|Я| = 5/2.

Область, в которой |Д| > 5/2, считается жидкой, а область, где |Д| < 5/2 — твёрдой. Внешняя граница области течения, описываемая формулами (2), выбирается так, чтобы жидкая зона целиком находилась между границами (1) и (2).

Пусть и, V — компоненты вектора скорости элементов матрицы по направлениям х, у. Введём функцию тока Ф и завихренность ш по формулам

_ 1 д Ф _ 1 д Ф _ 1 / д(иН) д ^Н) \ гН ду ' гН дх ' Н2 у ду дх у '

где Н = -у/(дг/дх)2 + (дг/ду)2 — коэффициент Ламэ.

Система уравнений для Ф, ш имеет вид [4]

дш

1

дЬ Н2

А (1(тушЛ + (1(туш)

дх \т дх ) ду \т ду

С,

д /1 дФ\ д /1 дФ\

дх \т дх ) ду \т ду )

Н 2ш,

где Ь — переменная времени,

С = — С = Н 2

д [ ду

дх дх

1 д ( 1 дФ\ 1 дНдФ

Н ду тН ду - тН3 дх дх

ду дУ

1д

1 д Ф\

1 дН дФ

Н ду \тН дх ) тН3 дх ду

д (ду / 1 д ( 1 дФ\ 1 дНдФ ду дх Н дх тН ду - тН3 ду дх

ду / 1 д / 1 дФ\ 1 дНдФ - ду Н дх тН дх тН3 ду ду

+

Ие / д /дФ

д дФ

(4)

(5)

2 \дх \дуШ) ду \дхШ Граничные условия системы следующие:

Ф = 0, дФ/ду = 0 при у = 0,

Ф = 0, ш = 0 при х = 0 и х = X,

Ф = -V (Ь)т2/2, ш = 0 при у = У. (7)

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

Ф = 0, ш = 0 V (Ь) = 0, при Ь = 0.

(8)

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

Руу = -Р, Рху = аш + ^п(ш),

(9)

где Р — гидродинамическое давление:

Г 1 д (туш)

р =----ах.

т ду

(10)

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

X

2

ат

X

2

сг1 = -- I Руут1 (х)ах, сХ2 = - I рт,ут1(х)ах

п

п

ах

(11)

где функции г^ж), ¿^(ж) определены формулами (1). Суммарная сила, действующая на тело, вычисляется по формуле

= + + Ао. (12)

Скорость тела определяется как

г

V (£) = ¿1 У (С^) + С^М + ^ (13)

о

а расстояние, пройденное телом, как

г

*!(*) = вJv (14)

о

Итак, ставится задача о решении системы уравнений (4), (5) с граничными условиями (7) и начальными условиями (8), где скорость тела V(£) определяется по формулам (9) — (13), а расстояние, пройденное телом, — по формуле (14).

4. Метод решения задачи

Уравнение (4) можно записать в виде

д N

+ ¿1 + М и = С, (15)

где

^ =-ндХ (ГдХ (г^), =-нК1ду (г^).

Пусть задана неравномерная сетка

0 = яо < Я1 < ... < -1 < жм = X, 0 = уо < У1 < ... < Ум-1 < Ум = У". Введём обозначения:

Хп —1/2 = (Хп-1 + Хп)/2, Л,ж„-1/2 = Хга - Хп-1, П =1,Х,

Ут-1/2 = (Ут-1 + Ут)/2, Лут-1/2 = Ут - Ут-1, т =1,М

Лхп = (Хп+1 - Хп-1)/2, п =1,х - 1, Лут = (Ут+1 - Ут-1 )/2, ТО =1,М - 1.

Будем считать, что

Лжга+1/2 - Лжга-1/2 = 0(ЛХп+1/2) , Лут+1/2 - Лут-1/2 = 0(ЛУт

Введём также равномерную сетку на полуоси £ > 0 с шагом т. Примем для любой функции /(я,У, £) обозначение

/ (яп, Ут, кт) /пт,

причем индексы п, т могут быть дробными либо отсутствовать, если f не зависит от х, у или эти индексы несущественны, а индекс к может отсутствовать, если f не зависит от Ь или этот индекс несуществен.

Аппроксимируем дифференциальные операторы Ь1, Ь2 разностными операторами Л1, Л2 следующим образом:

Л1шпт

Н 2 Ъ

Н пт ъхп

(тп+1 туп+1тшп+1т тптуптшпт)

тп+1/2тЪхп+1/2

(тпт уптшпт тп-1туп-1тшп-1т) тп 1 /2тЪхп-1 /2

Л2шпт = —

Н 2 ъ

Нптъут

( тпт+1 упт+1 шпт+1 тпт уптшпт )

тпт+1/2Ъут+1/2 (тптуптшпт тпт-1упт-1шпт—1)

тпт-1/2Ъут-1/2

Аппроксимируем (15) как

(шк+1 _ I к )

/I? I „,2~2 д к\к\ (шпт шпт) , / \к , д к\(, ,к+1 , ,к

(Е + 7 Т ЛкЛ2)- + 7(Л1 + Л2)(шпт - шпт) +

Т

+ (Лк + Лк)шкПт = Скпт, п =1,М - 1, т =1,М - 1, (16)

где Е — тождественный оператор, 7 € [0,1] — весовой коэффициент. Из уравнения (5) и первых двух условий (7) следует, что на границе тела приближённо выполняется условие Тома

2Фк+1

шк+1 = 2Фп1

тпо Н^о Ъу 1 / 2

Представим шпт в виде

ш ,к+1 = -к+1 + рк+1 ( 2Фп!1 _ ш к | (17)

шпт шпт Т 1 пт I тт2 Ь.2 по I > V ' 7

утпо НпоЪу1/2

где сеточная функция Р^т удовлетворяет задаче

(Е + 7тЛ2 ^ = 0, Ргако+1 = 1, Р^1 = 0, п = ^Г, а функция ^т1 — (16) и граничным условиям

¿к+1 = шпо, ¿Ц,1 = 0, п = 0,Ж ^ = ^И+т1 = 0, т = 6ТМ. (18)

Уравнение (16) можно записать в факторизованном виде

( +1 ) (Е + 1т Л к )(Е + 7т Л к )(шт-пт) = СкПт-

Т

-(Лк + Л к)шпт, п =1,Ж - 1, т =1,М - 1.

Отсюда следует, что функция О^т1 может быть вычислена по следующему алгоритму: 1) определяем невязку

Cm = Gm - (Л1 +Л2Km, П - 1, m = 1,M - 1;

2) для всех m = 1, M — 1 решаем системы уравнений относительно С

Com = 0,

k

nm

(E + 7ТЛ?)СПт = Cm, П = 1,N - 1, Cwm = 0;

3) при всех п = 1, N - 1 решаем системы уравнений относительно поправки

СПо = 0,

k .

ram"

(E + )Crim = С

m

1,M- 1,

С

nM

0,

полагаем

Com = 0, Ckm

0, m = 0, M;

4) для всех n = 0, N, m = 0, M вычисляем

k+1

+ rC

k

nm.

Подставляя в правую часть уравнения (5) выражение (17) и аппроксимируя левую часть (5) обычным образом, получим уравнение для функции тока (в силу (18) функцию

сП0 можно заменить на ССП+1)

h x

n+1m nm nm n—1m

rn+1/2mhxn+1/2 rn— 1/2mhxn—1/2

nm+1 nm nm nm— 1

+

h

ym

rnm+1/2hym+1/2 rnm—1/2hym—1/2 _

H2

.k+1

nm

, I - I P n

k+1 nm

n1

rn0 Hn0 hy1 / 2

k+1 nO

n = 1,N- 1, m = 1,M- 1,

которое можно записать в виде

^k+1 + bx ^k+1 + ay + by - c ^k+1+

г n— 1m 1 nm n+1m 1 nm nm— 1 1 nm nm+1 ^nm nm 1

+dnm+1^n+1 = -/nm1, n=t^n-i, m=1,m -1,

где

ax = nm

1

rn— 1/2m hxn—1/2 hx

=_1_

nm С

rnm—1/2 hym— 1 /2 hym

nm

1

by

nm

rn+1 / 2m hxn+1/2 hxn 1

rnm+1/2 hym+1 /2 h ym

k

1

1

x

a

с = 0х + Ьх + ау + ьу

пт пт пт пт пт 2 Н2 Рк+1

лук+1 _ _ пт пт г к+1 _ _ н 2 (д к+1 _ д к+1 р к+1)

^пт т н2 ъ2 , 3пт Нпт(шпт шпо Рпт ).

тпоНпоЪу1/2

К уравнению (19) нужно присоединить граничные условия (7):

Фп+1 = 0, Ф&1 = -Ук+1т2пМ/2, п = 6ТЖ, Фкт1 = Ф^ = 0, т = 0~М. (20)

Система уравнений (19), (20) отличается от системы (33), (34) работы [4] только отсутствием члена епт+1ФпМ^1 в левой части (19) и правыми частями в граничных условиях (20). Поэтому к ней можно полностью применить описанный в [4] способ решения, положив епт+1 = 0.

Замечание 1. Данный метод решения задачи пригоден только для случая, когда число Рейнольдса не велико. Если нелинейные конвективные члены в уравнении для вихря не являются малыми, так что число Рейнольдса велико, то нужно ввести модифицированную функцию вихря П = ш/т и решать задачу, как это сделано в работе [4].

5. Тесты

Все тестовые расчёты осуществлялись при числе Re = 0.

Для случая, когда форма тела является сферической, были произведены следующие тестовые расчёты:

1) при v = const численное стационарное решение исследовалось на сходимость к точному решению Стокса [5]. Установлен второй порядок сходимости вихря и функции тока к точным аналогам при измельчении расчётной сетки;

2) приближённо вычислялась функция |D| для решения Стокса. Для неё установлен второй порядок сходимости к точному аналогу при измельчении расчётной сетки;

3) вычислялись восемь членов (из которых шесть ненулевых) в выражении для функции G (см. формулу (6)) при vnm = 1 + xn + ym для решения Стокса. Для каждого из них установлен второй порядок сходимости к точным аналогам при измельчении расчётной сетки;

4) при vnm = 1 + x2n + ym исследовалось стационарное решение на сходимость "в себе" [4]. Установлен второй порядок сходимости вихря и функции тока при измельчении расчётной сетки;

5) при функции vnm, определяемой формулой (3), для 6 = 0.1, a =1 исследовалось стационарное решение на сходимость "в себе". Установлен второй порядок сходимости вихря и функции тока при измельчении расчётной сетки;

6) при vnm = 1 + Nb/(2|Dnm| + PR), Nb = 0.747, PR = 1.038 ■ 10-4 рассчитывалось стационарное решение, которое сравнивалось с решением, полученным в [6]. На рис. 1,1 приведены линии тока, полученные в настоящей работе, и одноименные линии тока, полученные для тех же данных в работе [6] (отличие знаков этих изолиний вызвано тем, что в [6] жидкость движется относительно капли в положительном, а в настоящей работе в отрицательном направлении относительно оси z — см. формулы (7), (20)). На рис. 1, II представлены изолинии функции |D|, полученные в настоящей работе и для тех же данных в работе [6]. Из этих рисунков видно, что результаты сравниваемых расчётов весьма близки, хотя полного совпадения нет, что можно объяснить меньшим числом разбиений в [6] (N = 16, M = 14) по сравнению с числом разбиений в настоящей работе (N = 80, M = 80);

Рис. 1. Расчётные линии тока при обтекании сферы жидкостью Бингама (I) и изолинии функции | (II), полученные в настоящей работе (а) и в [6] (б)

7) при v = const вычислялась скорость стационарного обтекания сферы и установлен второй порядок сходимости её величины к значению, описываемому точным решением Стокса.

Для эллиптической формы тела функции Ф(x,y,t), u(x,y,t), V(t) при обтекании эллипсоида жидкостью Бингама исследовались на сходимость "в себе" при t ^ то и для них был установлен искомый второй порядок сходимости при измельчении расчётной сетки.

Замечание 2. Для сферической формы тела в пунктах 1-6 скорость тела изменялась по заданному закону и стремилась к постоянному значению при Ь ^ то, в то время как в пункте 7 для сферической формы тела и в расчётах для эллиптической формы в каждый момент времени скорость тела определялась по формулам (9)-(13).

6. Пример расчёта

Рассмотрим класс эллипсоидов с большой полуосью длины то(Ь) и двумя малыми полуосями длины Ьто(Ь), где Ь < 1 — положительный числовой параметр. Величину то(Ь) подберём так, чтобы объём каждого эллипсоида из этого класса был равен объёму шара радиуса Тоо. Тогда

пл Тоо То(Ь)

Ь2/з

и формулы (1) в безразмерных переменных примут вид

r = b sin x, z = — cos x, x E [0,п]. Введём эллиптические координаты x, y :

r = p2 sin x sh (y + y0), z = —p2 cos x ch (y + y0), x E [0,n], y E [0,Y], (21)

где

P2

ch yo

b = \Л — p2 = th yo

Внешняя граница области течения будет описываться формулами (21) при у = У, х Е [0,п], части оси симметрии — этими же формулами при х = 0, х = п, у Е [0, У], граница тела — также этими формулами при у = 0, х Е [0, п]. Функции (21) конформно отображают прямоугольник 0 < х < п, 0 < у < У на область течения с коэффициентом

Ламэ _

Н = / (дт/дх)2 + (дт/ду)2 = р2 у'с^со2^.

Функция безразмерной вязкости матрицы при обтекании эллипсоида будет иметь вид (3) при

|Я| = /в2хх + Б2уу + Б* у + 2Б2ХЬ

где

DXX - " » ^

1

H

yy

1

H

D = — < - —

-ISiiii тт \ ^

D

д

dx

д dy

1

r2 H 2

1 д Ф rH dy

1 дФ

rH dx

дФ dr dy dx

+

xyi

1 dФ dH rH2 dx dy

1 дФ dH rH2 dy dx

д Ф dr

dx dy

D

xy

1 i d ' 1 5Ф" d ' 1 ЗФ"

2H [dy rH dy dx rH dx

1 dФ dH 1 дФ dH

+

rH2 dy dy rH2 dx dx

x E (0, n), y E (0, Y),

1

_ 1 М д 2Ф|

xx 2дт/дхду\Н2 дх2 f

d = _1А/ 1

yy H ду\ (дт/дх)Н дх2 J

}2v

D = дт/дх д f 1 д2Ф w = 2Н2 д^ \ (дт/дх)2 дх2

Dxy = 0, х = 0, х = п, у е (0, Y),

Dxx = 0, Dyy = 0, Dvv = 0, DXy = ш/2, у = 0, х е [0,п],

Dxx = 0, Dyy = 0, Dw = 0, Dxy = 0, у = Y, х е [0,п].

(Эти равенства верны не только для эллиптической, но и для любой достаточно гладкой формы тела.)

В таком случае формулы (11) примут вид

п п

2&2/3 > 26-1/3 Г

Gz1 =--Pyy sin х cos хвх, Gz2 =- Pxy sin2 хвх. (22)

п J п J

о о

Условием начала движения эллипсоида будет

A1 , ч

krit = AW' < 1 (23)

Данный критерий выполняется для нерегуляризованной модели жидкости Бингама, когда

1

V = a +

л/ад'

Действительно, учитывая, что на поверхности эллипсоида выполняются условия (9), (10), а завихренность ш в начальный отрезок времени мала, получим

Руу 0, Рху 1,

поскольку ш < 0. Подставляя эти значения в (22), имеем 0*1 = 0, 0*2 = -Ь-1/3. Требуя, чтобы в (13) доминировал последний член, получим (23).

При расчётах использовались следующие исходные данные:

и0 = 1.175 ■ 10-3 м2/с, р0 = 1260 кг/м3

(оба параметра соответствуют глицерину при 293 К [7]),

р1 = 7870 кг/м3

(железо [7]),

г00 = 0.01 м, к0 = 48.03 Па, = 0.9733 м/с, Ь0 = 0.2553 с, у1 = 3.917 ■ 10-4 м2/с,

У = 2.5, а = 3, Яе = 24.85, А1 = 0.3772, А2 = 2.161, В = 24.85, 5 = 0.0025, 7 = 0.999, N = 32, М = 32, е = 10-14 (последний параметр означает точность решения задачи для функции тока [4]).

Ниже приведены переменные, зависящие от Ь (здесь К — число итераций по времени):

Ь

Го (Ь),м кг11 т К

1

0.8

0.01 1.1604-10

-2

0.1745 10-4 4000

0.1880 10-4 4000

0.5 1.5874 • 10 0.2199 2•10-5 20 000

2

0.3 2.2314 • 10 0.2607 5•10-6 80 000

2

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

Число итераций при решении задачи для функции тока находилось в пределах от 1 до 150. Характерное значение максимума | составляло 0.25, что в сто раз больше величины 8.

На рис. 2, а показана расчётная сетка для Ь = 0.3 (соответствующая сетка в переменных х, у является равномерной), на рис. 2, б представлены форма расчётной области, в которой решалась задача (сплошная линия), и граница жидкой зоны (пунктир). На рис. 3 приведены зависимости скорости тела и расстояния, пройденного телом, от времени. Из рис. 3, а следует, что скорость стационарного обтекания эллипсоида зависит

Рис. 2. Расчётная сетка при Ь = 0.3 (а) и форма расчётной области (сплошная линия) и граница жидкой зоны в момент времени £ = 0.4 (пунктир) при Ь = 0.3 (б)

0.25г 0.20.150.10.050-

0 0.1 0.2 0.3 0.4

а

Рис. 3. Зависимость скорости тела (а) и расстояния, пройденного телом, (б) от времени. Сплошная линия — Ь =1, крупный пунктир — Ь = 0.8, мелкий пунктир — Ь = 0.5, штрихпунктир — Ь = 0.3

от Ь немонотонным образом: сначала с уменьшением Ь она снижается, а потом увеличивается. Кроме того, для Ь =1, 0.8 и 0.5 скорость тела в зависимости от времени является монотонно возрастающей функцией, в то время как для Ь = 0.3 она сначала возрастает, а потом монотонно убывает, стремясь к постоянному значению. Из рис. 3 также видно, что максимум скорости тела для всех вариантов близок к значению 0.02, т. е. примерно в 50 раз меньше скорости стационарного обтекания сферы ньютоновской жидкостью с вязкостью и0. Поэтому реальное число Рейнольдса Ие' = 0.02 Ие ~ 0.5. Установлен также следующий факт: для ньютоновской жидкости зависимость скорости стационарного обтекания эллипсоида от параметра Ь является монотонной (скорость снижается с уменьшением Ь). В общем случае здесь противопоставлены два фактора: с уменьшением Ь возрастает площадь поверхности эллипсоида. Это должно вести к снижению скорости обтекания, что и имеет место для случая ньютоновской жидкости, а также при достаточно больших Ь в случае жидкости Бингама. Вместе с тем при уменьшении Ь форма эллипсоида становится более обтекаемой, что может вызвать увеличение скорости обтекания, — последнее и наблюдается для жидкости Бингама при достаточно малом Ь.

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

Заключение

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

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

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

[1] Веретенцев В.А. Построение разностной сетки в области с криволинейными границами

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

М.: Изд-во Моск. гос. ун-та, 1989. С. 88-93.

[2] Пивоваров Ю.В. О построении ортогональной разностной сетки в криволинейном четырёхугольнике // Вычисл. технологии. 2003. Т. 8, № 5. С. 94-101.

[3] Malek J., Ruzicka M., Shelukhin V.V. Herschel — Bulkley fluids: Existence and regularity of steady flows // Math. Models and Methods in Appl. Sci. 2005. Vol. 15, No. 12. P. 1845-1861.

[4] Пивоваров Ю.В. Расчёт движения жидкости с переменной вязкостью в области с криволинейной границей // Вычисл. технологии. 2005. Т. 10, № 3. С. 87-107.

[5] Кочин Н.Е., Кибель И.А., Розе Н.В. Теоретическая гидромеханика. Часть II. М.: Гос. изд-во физ.-мат. лит-ры, 1963. 727 с.

[6] Beris A.N., Tsamopoulos J.A., Armstrong R.C., Brown R.A. Creeping motion of a sphere through a Bingham plastic // J. of Fluid Mechanics. 1985. Vol. 158, September. P. 219-244.

[7] Таблицы физических величин / Под ред. И.К. Кикоина. М.: Атомиздат, 1976. 1008 с.

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