CC BY
13
Вычислительные технологии
Том 22, № 3, 2017
Расчет движения двух твердых шаров в жидкости Бингама под действием силы тяжести
Ю. В. ПИВОВАРОВ
Институт гидродинамики им. М.А. Лаврентьева СО РАН, Новосибирск, Россия Контактный e-mail: pivov@hydro.nsc.ru
Рассматривается осесимметричное нестационарное движение двух одинаковых твердых шаров в жидкости Бингама под действием силы тяжести, направленной вдоль оси симметрии. В задаче используется подвижная ортогональная разностная сетка. Модель жидкости Бингама регуляризуется и сводится таким образом к описанию движения неньютоновской жидкости, у которой вязкость зависит определенным образом от второго инварианта тензора скоростей деформаций. Производятся расчеты для фиксированных свойств материалов жидкости и шаров в зависимости от радиуса шаров.
Ключевые слова: регуляризация, неньютоновская жидкость, конформное отображение, уравнения Навье — Стокса, завихренность, функция тока.
Введение
В природе, быту и технике существует широкий круг материалов, описываемых вязко-пластической моделью — моделью Бингама. Это такие материалы, как свежий бетон, глинистые почвы, нефтесодержащие материалы, сель, магма, коллоидные растворы, гели, порошкообразные смеси, смазочные материалы, металлы при обработке давлением, кровь в капиллярах, зубная паста, пищевые продукты (мучное тесто, мед, фруктовые сиропы и др.). Все эти материалы обладают общим свойством: ниже определенного предельного значения напряжений среда ведет себя как жесткое тело, выше этого предела — как несжимаемая вязкая жидкость.
Вязкопластические жидкости, как правило, неоднородны и часто содержат включения твердых частиц. В некоторых случаях эти включения желательно удалить, а в других — распределить в среде специальным способом для придания материалу нужных свойств. И в том и в другом случае необходимо знание особенностей движения и взаимодействия частиц в потоке. Например, изучение взаимодействия частиц важно для описания таких повсеместно встречающихся процессов, как осаждение твердых частиц и перенос газовых пузырьков.
Взаимодействию частиц в вязкой жидкости посвящено большое число работ (см., например, [1, 2] и цитируемую там литературу), в неньютоновской — гораздо меньше [3, 4], в вязкопластической среде — совсем мало. Ниже выделим некоторые работы о движении твердых тел в жидкости Бингама. Если движение погруженного тела является единственной причиной, вызывающей течение вязкопластической среды, то течение происходит только в ограниченной области вокруг тела, а остальной объем остается неподвижен. Эта особенность отличает движение включений в вязкопластических
© ИВТ СО РАН, 2017
средах от движения в других жидкостях. Определение области течения вокруг тела является одной из основных целей исследований в данной области. В [5] сделаны расчеты движения шара при нулевом числе Рейнольдса. Было получено, что движение шара осуществляется только при числе Бингама, меньшем определенной величины. При пределе текучести, близком к соответствующему значению, около шара образуется пластический пограничный слой. В [6] повторены некоторые результаты [5] и сделаны расчеты для эллиптической формы тела и числа Рейнольдса порядка единицы. Был рассмотрен класс осесимметричных эллипсоидов, имеющих одинаковый объем, но разное отношение полуосей. Показано, что, во-первых, зависимость скорости стационарного обтекания от отношения полуосей эллипсоида носит немонотонный характер, а во-вторых, что при достаточно вытянутой форме эллипсоида его скорость в зависимости от времени сначала увеличивается, а потом уменьшается, стремясь к постоянному пределу. В [7] экспериментально изучено медленное движение одного и двух шаров. Определена зависимость коэффициента сопротивления от шероховатости поверхности шаров и расстояния между ними. В [8] сделаны расчеты, описывающие взаимодействие двух шаров, движущихся вдоль линии, соединяющей их центры. Показано, что для безразмерного расстояния между центрами шаров, большего шести, жидкие зоны вокруг шаров идентичны их аналогу вокруг одного шара.
Отметим, что в [7, 8] рассматривалось только стационарное так называемое ползущее течение вокруг шаров, т. е. течение при отсутствии конвективных членов в уравнении импульса. В этом случае оно является симметричным относительно плоскости, проходящей через середину отрезка, соединяющего центры шаров, ортогонально ему. Шары при этом движутся с одинаковыми скоростями. В экспериментах по осаждению двух одинаковых частиц под действием силы тяжести, как правило, скорость первой заметно ниже, чем второй, что делает течение нестационарным. Одной из причин этого эффекта может быть инерционное взаимодействие, моделирование которого требует учета конвективных членов в уравнении импульса.
В настоящей работе рассматривается осесимметричное нестационарное движение двух одинаковых твердых шаров в жидкости Бингама под действием силы тяжести, направленной вдоль оси симметрии с учетом конвективных членов в уравнении импульса, наличие которых обусловливает движение шаров друг относительно друга. Источником движения является сила тяжести. Как отмечалось выше, в каждый момент времени вне некоторой поверхности, охватывающей шары, поле скоростей частиц жидкости тождественно равно нулю. В цилиндрической системе координат r,z,ip с началом в середине отрезка, соединяющего центры шаров, где г — полярный радиус, ip — полярный угол, от которого искомые функции не зависят, удобно решать задачу в полуплоскости г > 0, ip = const. При определенном выборе формы охватывающей шары поверхности область течения в переменных г, z является криволинейным четырехугольником. Она будет двигаться вместе с шарами и одновременно деформироваться, так как шары движутся друг относительно друга.
Задачу предлагается решать конечно-разностным методом с использованием подвижной ортогональной разностной сетки. При построении последней предварительно насчитывается 15 так называемых базовых сеток, отвечающих дискретному набору расстояний между центрами шаров. Для каждого текущего значения расстояния между центрами шаров граничные узлы сетки строятся методом линейной интерполяции с граничных узлов базовых сеток, соответствующих двум соседним значениям из дискретного набора расстояний между центрами шаров, а внутренние узлы сетки нахо-
дятся из решения уравнений Лапласа для функций гсога/(х, у, I), zconf (х, у, I), где х,у — переменные, принадлежащие прямоугольной области, которую функции rconf, zconf конформно отображают на область течения в переменных г,х, Ь — переменная времени.
При решении задачи модель жидкости Бингама регуляризуется и сводится таким образом к модели неньютоновской жидкости, у которой вязкость зависит определенным образом от второго инварианта тензора скоростей деформаций. В [9] доказано, что при параметре регуляризации, стремящемся к нулю, решение регуляризованной задачи сходится к решению исходной задачи для жидкости Бингама. Предлагаемый подход к решению задачи оригинален и более прост (за исключением метода построения расчетной сетки), чем метод конечных элементов, используемый в [8], поскольку в первом после расщепления по направлениям приходится решать системы линейных алгебраических уравнений с трех- и пятидиагональными матрицами, для которых существует простой прямой метод решения — метод прогонки, тогда как в последнем обращаемая матрица хотя и сильно разрежена, но имеет более сложную структуру. Кроме того, метод, используемый в [8], является более громоздким, более сложен в программировании, требует больше памяти и вычислительных затрат.
1. Построение подвижной ортогональной разностной сетки
Пусть г1 — расстояние между центрами шаров, за единицу измерения которого принят радиус шаров г0. В каждый момент времени движения шаров ортогональная разностная сетка в области течения зависит только от г1. Будем считать, что г1 Е [2.01, 9]. На рис. 1 показана часть И области течения, лежащая в полуплоскости г > 0. Здесь = Т\/2 — координата г центра правого шара, Г! — граница правого шара, Г2 — часть оси симметрии перед правым шаром, Г3 — проекция плоскости симметрии между шарами на плоскость рисунка, Г4 — кривая, частично совпадающая с осью симметрии Ох и подходящая под прямым углом к отрезку Г3. Она состоит из трех участков: прямолинейного участка длиной I = 4, 1/4 окружности радиусом г3 = 8 с центром в точке г = г3, г = г1 /2 + 1 + I и 1/4 окружности радиусом г2 = г1/2 + 1 + I + г3 с центром в точке г = г3, г = 0 (в качестве масштаба длины используется радиус шара г0). Вектор силы тяжести направлен справа налево.
г
Рис. 1. Схема задачи и система координат
Первым этапом решения задачи о движении шаров в жидкости Бингама является построение сеток для г1 е (2.01, 2.02, 2.04, 2.07, 2.1, 2.2, 2.4, 2.6, 3, 4, 5,6, 7, 8, 9}. Эти сетки строятся до начала расчета движения шаров и используются потом для построения сеток с произвольным г1 из заданного промежутка. Поэтому будем называть их базовыми. Способ построения этих сеток подробно описан в работе [10].
Обозначим через П прямоугольную область 0 < х < 1, 0 < у < У/2, которая конформно отображается на часть И области течения, находящуюся в полуплоскости г > 0, с помощью функций rconf (х,у,Ь), (х,у,Ь). Причем отрезок 0 < х < 1, у = 0 переходит в границу Г!, отрезок х = 0, 0 < у < У/2 — в границу Г2, отрезок 0 < х < 1, у = У/2 — в границу Г3, отрезок х =1, 0 < у < У/2 — в границу Г4. Пусть 0 = х0 < Х\ < ... < Хм-1 < хм = 1, 0 = у0 < у\ < ... < ум/2-1 < Ум/2 = У/2 — одномерные сетки в прямоугольнике П, причем М, М делятся на 4. Образом этих сеток в области И будут ортогональные сетки = гСоп/(хп,ут,Ьк), = (хп,ут,Ьк), п = 0, И, т =
0, М/2, где к — номер итерации по времени.
Опишем свойства, которым должна удовлетворять одномерная сетка хп, п = 0,Ы. Обозначим
кп-1/2
Хп Хп—1,
П
и потребуем, чтобы
к1/2
к
11 шт
< 1/N,
-1/2
к
12 шт
< 1/N,
кп+1/2/кп-1/2
еопя1 > 1,
п = 1, N/2 — 1, кМ/2-1/2 = кх
N/2+1/2,
, К-1/2/кп+1/2 = еОП^ > 1
п = Ы/2 + 1,Ы — 1.
Этими условиями сетка хп, п = 0,Ы, определяется однозначно. Алгоритм построения такой сетки подробно изложен в [10]. В качестве исходных данных необходимо задать число разбиений М, а также первый к11ш;п и последний к12ш{п шаги сетки. Задавая число разбиений М', первый к21ш\п и последний к22тт шаги сетки у'т, построим эту сетку так же, как строили выше сетку хп. Далее, полагая М' = М/2, сетку ут определяем по формуле
Ут = у'тХ/2, т = 0,Ж
Введем в рассмотрение геометрические параметры базовых сеток в зависимости от г1: среднее вас и минимальное ва ш;п значения синуса угла между координатными линиями:
вт а[г
гхГу гуТх
(^Х + ГХ )(гу + Гу)
п = 0,^- 1, т = 0,М'~ 1,
где
I _ %п+1,т ^п,т | _ %п,т+1 %п,т
%хуп,т ~ ~ , 1п,т
Хп+1 ХГ:
Ут+1 Ут
формулы для гх1п,т, Ту |п,т выглядят аналогично, далее
М-1 м'-1
^ ^ вЬ а[п,т к
хп+1/2кут+1/2
п=0 т=0
/
м1
М'-1
^^ кхп+1/2 У ] кут+1/2
п=0
т=0
hxn+1/2 %п+1 hym+1/2 Ут+1 Уп
Sa min = min{sin о\ : п = 0,N — 1, т = 0, М' — 1}.
Положим
SL = (1 — Sac) • 100 %, S*a max = (1 — Sa min) • 100 %
— отличия от единицы среднего и минимального значений sin a.
В табл. 1 приведены значения Y, значения h12min, умноженные на N, значения S*c, S*max и значения относительных среднеквадратических отклонений от нуля ошибок в выполнении условий Коши — Римана а1, о2. Видно, что наихудшим качеством обладает первая сетка. Это объясняется очень малым расстоянием между перифериями шаров для этой сетки. Качество остальных сеток удовлетворительное.
Перейдем к описанию процесса построения сеток при произвольном г1 Е [2.01, 9]. Пусть г\а и Г\в — соседние значения г1 из дискретного набора, а текущее значение г1 таково, что
Г\А < П < пв.
Тогда любую числовую характеристику X сетки из описанного ниже множества характеристик определим с помощью формулы линейной интерполяции
х = (п — Г1А)ХВ + (пв — п)ХА Г\в — Па '
Здесь Ха, Хв — характеристика X базовых сеток с г\ = г\а и г\ = Г\В соответственно. В качестве X будем принимать следующие характеристики:
а) значение Y;
б) значение hi2min;
в) значения rn,o, znfl, rn,M, zn,M, п = 0,N; Го,т, Zo,m, rN,m, zN,m, т =1,М — 1. Значения h11min = 1/N, h21min = 1/М', h22min = 0.5/М' являются одинаковыми для
всех сеток.
Таблица 1. Зависимость параметров сетки от п
к Г1 Y $ас, % с * % °a max, % % О 2, %
1 2.01 0.069 720 0.03 125 0.334 5.73 23.5 9.37
2 2.02 0.088 852 0.03 125 0.061 5.07 15.3 7.97
3 2.04 0.133315 0.0 625 0.091 3.96 16.6 9.32
4 2.07 0.174 386 0.0 625 0.115 3.52 15.1 10.1
5 2.10 0.208 080 0.0 625 0.568 4.97 16.2 12.6
6 2.20 0.282 759 0.0 625 0.177 2.94 11.5 11.1
7 2.40 0.395 726 0.125 0.156 3.21 9.56 9.50
8 2.60 0.482 370 0.125 0.181 2.98 9.03 9.57
9 3.00 0.615 503 0.125 0.210 2.75 8.44 9.40
10 4.00 0.847017 0.25 0.201 2.96 7.14 8.11
11 5.00 1.012 498 0.25 0.217 2.73 6.99 7.95
12 6.00 1.144 376 0.25 0.224 2.78 6.93 7.82
13 7.00 1.255 358 0.25 0.229 2.78 6.90 7.73
14 8.00 1.351923 0.25 0.233 2.73 6.90 7.65
15 9.00 1.437629 0.5 0.211 3.17 6.45 7.22
Рис. 2. Прямоугольная сетка в переменных х, у (а) и ортогональная разностная сетка в области Их в переменных г, г при гх = 3 (б)
С помощью формул
ут :=2ум/2 - Ум-т, т = М/2 + 1,М;
гп,т := гп,м-ш, 2п,т := -гпМ-т, п = 0,М, т = М/2 + 1,М,
к области И присоединяется ее образ в полуплоскости г < 0. Полученную область обозначим .
На рис. 2 показаны базовая сетка в области и ее прообраз в переменных х,у при гх = 3. Видно, что она сгущается вблизи поверхности шаров и разреживается при удалении от них.
Суть метода построения ортогональной сетки для текущей итерации по времени состоит в следующем. С помощью известной ортогональной сетки на предыдущей итерации и неортогональных сеток на предыдущей и текущей итерациях, строящихся методом трансфинитной интерполяции с границы [11, 12], мало отличающихся друг от друга в силу малости шага по времени, конструируется неортогональная сетка, близкая к ортогональной сетке на предыдущей итерации. Она удовлетворяет всем необходимым граничным условиям, а значит, является хорошим начальным приближением для построения ортогональной сетки на текущей итерации. После построения начального приближения решаются разностные аналоги уравнений Лапласа для функций Гсоп/(%п,Ут^к), ¿соп/(%п,Ут^к), которые и определяют искомую ортогональную сетку на текущей итерации по времени.
2. Параметры задачи и искомые функции
Размерными параметрами задачи являются: ускорение свободного падения д = 9.81 м/с2, кинематическая вязкость материала матрицы ^ = 1.175в-3 м2/с (глицерин), плотность материала матрицы ро = 1260 кг/м3 (глицерин [13]), плотность материала шаров рх = 7870 кг/м3 (железо), радиус шаров гоо = 0.01 м, характерная скорость шаров (формула Стокса [14])
Уоо = 2дг1о ——— = 0.9733 м/с, 3 щр0
предел текучести материала матрицы к0 = 48.033 Па (выбирается так, чтобы описываемый ниже безразмерный параметр а0 был равен трем), характерная вязкость, определяемая пределом текучести материала матрицы,
^о = ^ = 3.917е-4 м2/с, РоУоо
масса шара
4 3
m10 = -пr00pi = 3.297е-2 кг, 3
характерный промежуток времени
t00 = wV00 = 1.0274е- 2 с.
Безразмерные параметры задачи: обратное число Бингама
«0 = щ/^10 = 3,
число Рейнольдса
Re0 = V00W^10 = 24.85, отношение характерной силы к произведению массы шара на его ускорение
Ли, = t00n\°°k0 = 1.5180е-2, m10 V00
отношение разности силы тяжести и силы плавучести к произведению массы шара на его ускорение
л t00(P 10 - 900)9 о ЙПо о А20 =--= 8.698е-2.
Pl0 V00
Описанные выше размерные и безразмерные параметры отвечают значению радиуса шаров г00 = 0.01 м. Их значения для произвольного радиуса шаров г0 определяются по формулам (в обозначении параметров снят нижний правый индекс 0):
V0 = V00C2, vi = uw/(, mi = mw(3, h = WC, a = «0C, Re = Re0C4, Ai = Aw/(4, A = М0/С3, < = VГ00.
Искомыми функциями являются: компоненты скорости u, v в направлениях х, у; функция тока Ф и завихренность ш, связанные с компонентами скорости и коэффициентом Ламе Н следующим образом:
1д Ф 1д Ф 1 (д (иН) д (vH)
и = —— -—, V =--— -—, ш
(д (иН) д (уН ) \ г Н ду ' гН дх ' Н2 \ ду дх ) '
Н = (дг/дх)2 + (дг/ду )2;
компоненты тензора скоростей деформаций И (они отличаются от аналогичных компонент в работе [6, с. 72-73] только тем, что граница х = ж заменена на границу х = 1, а значения компоненты Иху = 0 при у = V — на Иху = ш/2); модуль (корень из второго инварианта) тензора скоростей деформаций
|И| = + , х е [0,1], уе [0,V];
функция вязкости материала матрицы V(|^|); гидродинамическое давление рг; касательные Р^у и нормальные Руу напряжения на поверхности шаров; составляющие сил, действующих на шары, обусловленные нормальными Сгг1 и касательными Сгх2 напряжениями; суммарные силы С\, действующие на шары; скорости шаров Vг; расстояния гг, пройденные шарами, причем г = Ь соответствует левому, г = Р — правому шару.
3. Постановка задачи
Регуляризованная модель жидкости Бингама является моделью неньютоновской жидкости, у которой вязкость зависит от модуля тензора скоростей деформаций следующим образом:
Г а + 1/(^ВД), |Я| > 38/2, и = I f (т, 8/2 38/2,
( а + 1/(^28), ^ <8/2.
Здесь 8 — параметр регуляризации (малое число); /(|D|) — полином пятой степени, осуществляющий "склейку" значений функции и(|^|), ее первых и вторых производных в точках ^ = 8/2 и ^ = 38/2 [6].
Жидкая зона занимает область, описываемую соотношением
^ > 8/2.
Область, для которой |D| < 8/2, считается твердой.
Система уравнений для функции тока Ф и завихренности ш имеет вид
(1)
дш 1 ~dt + Н2
ЗФ дш ду дх
ЗФ дш дх ду
1
Re
где t — переменная времени,
G1
2
dv ду
Re Я 2
f 1 (дф \
\ Н ду \ гН дх )
(-(
\дх \
А (1дФ \ А (1 н2
дх \г дх ) ду \г ду )
G = G1 + G2 + G3 f д \dv ( 1 д / 1 дФ\ \дх дх \ Н ду \ гН ду )
1А (гииЛ + А А (гиш)
г дх ) ду\гду
G., (2) (3)
+
д_ дх
1
dv ( 1 д дх\ Н ду
дНд Ф
гН3 дх ду
д_ ду
1 д Ф гН ду
dv
дх V Н
1 дНдФ
гН3 дх дх i 1 д ( 1 д Ф \Ндх V
)
)
1 дНд Ф
гН ду J гН3 ду дх
)
dv f 1 д / 1 дФ\ 1 дНдФ ду \Н дх \гН дх ) гН3 ду ду
G2
дш дх дх dt
G3 = -
2шд_Н ~E~dt
+ н 2
д_ дх
(
дш ду ду dt'
1 дНдФ гН dt dx
) + ду('>
1 dHd Ф
гН dt dy
Здесь член G1 обусловлен непостоянством вязкости v, член G2 — зависимостью от времени сеток хп, ут, член G3 — неперестановочностью операторов rot и d/dt для вектора скорости (u,v).
Заметим, что конвективные члены в уравнении (2) записаны в недивергентной форме. Это связано с тем, что они будут аппроксимироваться разностями против потока второго порядка точности, которые не являются консервативными даже для дивергентной формы записи этих членов. Это сделано для того, чтобы уменьшить влияние схемной вязкости в алгоритме.
Граничные условия для функции тока и завихренности следующие:
„ т Ур -Уь 2 дФ Ур -Уь дг у = 0 : Ф =--г , -— =--г——;
4 ду 2 ду
х = 0 : Ф = 0, ш = 0; (4)
т Ур -Уь 2 дФ Ур -Уь дг у = ¥ :Ф = ^^Г, = —^гдТу;
Ур + УЬ 2
х = 1: Ф=---г2, ш = 0
4
(система координат движется со скоростью ( Ур + Уь)/22). Здесь первая пара условий — условия прилипания на границе правого шара, вторая — условия симметрии на оси симметрии между шарами, третья пара — условия прилипания на границе левого шара, последняя пара — условия симметрии на прямолинейном участке границы, т. е. при г = 0, и условия твердотельного движения со скоростью ( Ур + Уь)/2 на криволинейном участке внешней границы, т. е. при г > 0.
Начальные условия
Ь = 0 : Ф = 0, ш = 0, Уь = 0, Ур = 0. (5)
Функции гидродинамического давления на поверхности шаров определяются по формуле
х
*}■ (6)
0
Выражения для нормальных и касательных напряжений на границах шаров имеют вид
Ргуу = -Рг, Ргху = "ш, ге[Ь,Р}. (7)
Составляющие силы, обусловленные нормальными и касательными напряжениями в гидродинамическом течении вокруг шаров, равны
1 1
/р дт р [' р дт
РуудХГ(1х, = -2] Рху~дуГС^Х,
оо
1 1
/Е/ дг ^ [' дт
Руу~оХГ(Сх, = 2у Рху~дуГ(Сх. (8)
оо
Суммарная сила, действующая на -й шар, вычисляется по формуле
с; = 0^ + с;2 + А2/Ах , ге{Ь,Р}. (9)
Скорость г-го шара определяется как
V*(¿) = (£)£, г Е [Ь, Р}, (10)
о
а расстояние, пройденное г-м шаром, как
х\г) = I Г (1)£, г Е [Ь,Р}. о
Расстояние между центрами шаров
П(1) = по - (гр(*) — ^(¿)),
где Г10 — начальное значение г1.
Итак, ставится задача о решении системы уравнений (2), (3) с граничными (4) и начальными (5) условиями, где скорости Уь (¿) и Ур (¿) определяются по формулам (6)-(10).
4. Метод решения задачи
Уравнение (2) можно записать в виде
(д/дг + Ь + Ь2 + Ьз + ЬА)ш
где
с,
Ь1ш = —
1 д
Н2И,е дх \ г дх
(1д ( ) I -— (гиш)
,
Ь2ш = —
1 д
Н2И,е ду \г ду
(1д ( ) I -— (гиш)
)
Ьш_±_№дш Ьш- д!дш
Н2 ду дх' Н2 дх ду
Пусть к — номер итерации по времени, Тк переменной Ь на к-й итерации,
к-й шаг по времени, ^ — значение
* и = Ё
Т<п.
п=1
Примем для любой функции ¡(х,у, ¿) обозначение
!(хп, Ут, ^к) I'п,т,
причем индексы п,т могут быть дробными либо отсутствовать, если £ не зависит от х, у или эти индексы не существенны, а индекс к может отсутствовать, если £ не зависит от или этот индекс не существен.
Аппроксимируем дифференциальные операторы Ь1,Ь2,Ь3,Ь4 разностными операторами Л1, Л2, Л3, Л4. Первые два оператора аппроксимируются так же, как в [6, с. 68]. Операторы Л3, Л4 аппроксимируются следующим образом:
Л3шп,1
1
Н 2 1ап
1 п,т
['с1пи,п,тшп—2,т ^>пип,тшп-1,т+
где
+ (Сп ип,т Сп ипШп,т +с1пип,тШп+1,т еп^п,тШп+2,т] ,
Л4Шп,т ^ \Р'тУ^п,т^п,т—2 ^тУ^п,т^п,т— 1 +
1 п,т
+ {Ст ^п,т - СУп шп,т +(ЩГУ п,тШ п,т+1 - еупУп,тШп,т+2^ ,
х _ К-1/2 ,х _ 2^п~ 1 ]Х _ 2^п+1 рх _ ^1/2
,х __п~1/2 ^х - ---п-1 1Х - ---п+1 рх
п _ о^х Ъх ' п _ ъх ъх , п _ ^х ^Х ' п 0ъх ЪХ
п+
2Ъх Ьх п Ьх Ьх п Ьх Ьх п 2Ъх Ьх
2Пп-3/2 п— 1 Ьп-1/2 п— 3/2 Ьп+1/2Ьп+3/2 2Ьп+3/2Ьп+1
°п _ К - 0>п, сп _ С - е'п, п _2,М - 2, т _1,М - 1
и^. /т, _
Фп,т+1 ^п,т-1 ~т~ — Vп,т |^п,т| ^-1-+ Vп,т + |^п,т|
~ Т 11 1 ^ п гп 1 ^ п гп
2Ьт ' ,т 2 ' п,т 2
Фп+1,т Фп— 1,т
Упт _--п+1,т -1,т, п _1,М - 1, т _1,М - 1,
2Ьх
К-1/2 _ хп - Хп-1, п _1,^, кХп _ (Хп+1 -Хп-1)/2, п - 1,
Кг-1/2 _ Ут - Ут-1, т _1,М, НУт _(Ут+1 - Ут-1)/2, т _ 1, М - 1.
Формулы для У~т, У,+т, аут, ьут, (Щп, е~т, сУт, сУ+ выглядят аналогично.
Таким образом вязкие и конвективные члены аппроксимируются с вторым порядком, причем для конвективных членов вычисляются разности против потока. Это делается с целью уменьшения влияния схемной вязкости в алгоритме. В приграничных узлах сетки разности в сторону границы берутся с первым порядком аппроксимации. Вследствие этого и весь алгоритм теоретически будет иметь первый порядок. Обозначив
,т (Л1 + Л3 К ,т Л^Шп ,т (Л2 + Л4)и п,т аппроксимируем (2) как ^.к+1 _
(Е + 72 г2 Л5 Лд ) -^ + 7(Л^ + Лд )(ип+т - <т) + (Л5 + Лд )<„ _ С^, (11)
Тк , , , ,
где Е — тождественный оператор, 7 Е [0,1] — весовой коэффициент.
Из уравнения (3) и условий (4) для у _ 0 и у _ V следует, что на границах шаров приближенно выполняются условия Тома:
Ь+1 _ 2Фп,1 + рР ь+1 _ 2Фп,М—1 +
ШпА _ г Ч2 Ьу2 + Рп , Шп,М _ г Ч2 Ьу2 + Рп ,
I п,011п,0Ь1/2 'п,М^п,МЬМ —1/2
Р_УР -уь ( <1
с
Рп о~, и2 I 7 у2 + (Гп,1 гп,о)Гхх
2гп,сЯп,е \
уР -Vе I г2п
п,0 ,
Т^Ь у - у I ' п,М—1 , ч
Рп _ - -^^ 7^2--(Гп,М - Гп,М—1)ГХ
2Гп,М^п,М \ЬУМ—1/2
п, М
Представим шп+ в виде
,к+1 = шк+1 I рк+1 п, т п, т
2ФН1
п,1
,,2 , у2 г 1 п — шп,0 Тп.оП^ /О
р к+1
+ К- ш
0Нп,0^1/2
+
к+1
п, т
2Фк+
1
п,М-1 , рЬ шк+1
г рп — шп М
_гп,МНп,М Щ2-1/2
(12)
где сеточные функции Рк+, ОЛУТ есть решения следующих краевых задач:
(Е + 1тЛ к )Рк+1 = 0, Рк+о1 = 1, Рпк+ = 0, (Е + 7гЛб )Яп+т = 0, яп+о1 = 0, яп+м = 1,
а сеточная функция шп+ удовлетворяет уравнению (11) и граничным условиям
шп++1 = шп,о, шп+М = шкпМ, п = 1,Ж — 1, шк+т1 = ш^+т = 0, т = б;ж
Уравнение (11) можно записать в факторизованном виде
ш к+1 ш к
(Е + 7 тк Л к )(Е + 7 гк Л к) -^ = С^т — (Л5 + Лк
к
п = 1,# — 1, т =1,М — 1.
ш к+1 ш к п, т — п, т
п, т п, т
Причем граничными условиями для поправки С,пт =-являются однородные
к
условия первого рода.
Заметим, что на втором и третьем шагах алгоритма решаются линейные алгебраические системы уравнений с пятидиагональными матрицами. Шаг по времени Тк на к-м временном слое выбирается так, чтобы они имели диагональное преобладание с запасом не менее 25 %.
Подставляя выражение (12) в правую часть второго уравнения (2) и аппроксимируя левую его часть обычным образом, получим разностное уравнение для функции тока:
^птФп- 1,т + ^птФп+ 1,т + ^птФп,т-1 + ЬптФп,т-\-1 Сп,т!п,т+
+<П 1Фп,1 + еупкт 1Фп, м-1 = — и,т, п =1,М — 1, т =1,М — 1, (13)
где ^пт Кт апт <т — те же, что и в [6, с. 69-70К а
2 Н 2 С)к+1
у 2 Нп,тя п,т
^п,т
г п,МНп,М Щ^-1/2
г = _ Н2 Гшк+1 + рк+1 ГрР _ шк+1] + Пк+1 ГрЬ _ шк+111 }п,т Н п,т \шп,т + Р п,т п шп,о ] + яп,т п шп,М\ ] .
К уравнению (13) нужно присоединить граничные условия (4):
т Vр — VЬ 2 т Vр — VЬ 2 —_
Фп,о =--4-rn,о, фп,М = -4-Гп ,М, п =1^ — ^
V р + V Ь _
Фо,т = 0, Фм,т =-4-г%т, т = 0,М. (14)
Метод решения системы уравнений (13), (14) подробно изложен в работе [15].
5. Результаты численных экспериментов
Используемый численный алгоритм разработан на основе его аналога для задачи о движении одного шара в жидкости Бингама. Для этого аналога в [6] проведены многочисленные тесты. В том числе осуществлено сравнение с точным решением и результатами расчета других авторов [5]. Для нового элемента — метода пятидиагональной прогонки в уравнении для завихренности — произведены соответствующие тестовые расчеты. Точность решения задачи для функций rcoпf (х, у), гсоп/(х, у) при построении конформного отображения прямоугольника на область течения контролировалась в процессе вычислений.
При расчетах использовались следующие исходные данные: числа разбиений по оси хЫ _ 24 и оси уМ _ 48, число итераций при построении подвижной ортогональной разностной сетки Кс _ 2. При этом точность в решении уравнений Лапласа для функций гсоп/ (х, у), гсоп/(х, у) составляет 10—12 - 10—3. Значение параметра регуляризации 8 _ 0.02, значение весового коэффициента 7 _ 0.999, точность решения уравнения для функции тока
£ Ф
\
;Фпя ^ "ял+///2Я 2Й<10—9
где Фп — функция невязки; / — функция правой части в уравнении (13); безразмерный начальный шаг по времени 0 _ 10—3. Остальные параметры подробно описаны в разд. 2.
Сначала были произведены расчеты при радиусе шаров г0 _ 0.005 м и различных начальных значениях безразмерного расстояния между центрами шаров г10 Е [2.5, 9]. При этом
с _ 0.5, У0 _ 0.2433 м/с, ^ _ 7.833е-4 м2/с, т1 _ 4.121е-3 кг, ¿0 _ 2.055е-2 с, а _ 1.5, И,е _ 1.5531, Д _ 0.2429, А2 _ 0.6958.
Число итераций по времени составляло Ктах _ 20 000.
На рис. 3 показаны зависимости скорости левого шара от времени для г10 Е {2.5, 3, 4, 5, 6, 7, 8, 9}. Для всех вариантов суммарные силы СЬ, Ср, действующие на шары, с течением времени экспоненциально стремились к нулю, так что на заключительном этапе расчета они составляли порядка 10—10-10—8. Так как формальное число Рейнольдса И,е ~ 1.5, а характерное значение безразмерной скорости Vе ~ 0.02, то фактическое значение числа Рейнольдса И,е / ~ 0.03, т.е. является малым. Это объясняет тот факт, что расстояние между шарами меняется, поскольку скорость правого шара немного больше, чем левого, но это изменение происходит очень медленно. Ясно, что поскольку суммарные силы, действующие на шары при номере итераций по времени к > Ктах, отсутствуют, шары будут сближаться с постоянной скоростью, пока не столкнутся.
Рассчитаем время , за которое расстояние между центрами шаров уменьшается на величину их радиуса г0, и время от начала движения до соударения шаров ¿*:
** _ У0(УР°- УЬ), Г _ П(Г 10 - 2).
V, м/с
0.004 -
0.002 -
0
О 0.1 0.2 0.3 с
Рис. 3. Зависимость скорости движения левого шара от времени при го = 0.005 м и различных значениях г\0 € {2.5, 3, 4, 5, 6, 7, 8, 9} — кривые 1-8 соответственно
В табл. 2 приведены значения Vр, Vь, АУ = Vо(Vр — Vь), ¿* в зависимости от начального расстояния между центрами шаров г 10. Видно, что время сближения шаров на величину г0 составляет от 15 до 30 мин, а время от начала движения до соударения шаров — от 15 мин до 1 ч 45 мин.
Необходимо отметить следующее: если рассчитать форму жидких зон вокруг обоих шаров для описанных выше данных, то увидим, что они не только не пересекаются, но даже не полностью покрывают поверхности шаров. Это противоречие объясняется дефектом регуляризованной модели. (Если, например, жидкую зону определить неравенством |Д| > 5/4, то она будет охватывать поверхности шаров.) Для данных, описанных ниже, жидкая зона, определяемая неравенством (1), охватывает оба шара.
Вторая серия расчетов была произведена для безразмерного расстояния между центрами шаров г 10 = 3 и различных значений радиуса шаров г0.
В табл. 3 приведены характеристики течения в зависимости от радиуса шаров о. При этом величины, у которых размерность не указана, являются безразмерными. Шаг по времени т брался равным минимуму из двух значений: т0 иг* • 0.75, где т0 = 10-3, а т* — значение т, обеспечивающее на втором и третьем этапах решения задачи для ш (см. разд. 4) диагональное преобладание в матрицах Е+и Е+^тЛ^, Ктах — общее число итераций по времени, Vь(Ктах), Vр (Хтах), т(Ктах) — значения Vь, Vр, г на заключительном этапе расчетов ( к = Хтах). Фактическое значение числа Рейнольдса
Таблица 2. Зависимость параметров задачи от г\0 при г0 = 0.005 м
п о Vр • 102 Vь • 102 ДУ • 106, м/с **, с Г, с
2.5 2.06 956 2.06 839 2.847 11756 878
3 2.01454 2.01301 3.722 1 343 1343
4 1.85 521 1.85 323 4.817 1 038 2 076
5 1.79 279 1.79 083 4.769 1 049 3145
6 1.69 790 1.69 563 5.523 905 3 621
7 1.66 036 1.65 802 5.693 878 4 391
8 1.63 761 1.63 566 4.744 1 054 6 323
9 1.59 351 1.59126 5.474 913 6 394
12 3 4
5 6 7 8
Таблица 3. Зависимость параметров задачи от г0 при по = 3
Го, м 0.006 0.007 0.0085 0.01 0.02
с 0.6 0.7 0.85 1 2
^тах • Ю 3 200 2 800 2 400 2 000 1000
^(^тах) • 102 7.171 6.472 6.367 6.671 12.516
(^тах) • 102 7.213 6.499 6.408 6.709 12.592
Г(Ктах) • 105 10.60 14.12 9.79 2.07 3.63
Ус, м/с 0.3 504 0.4 769 0.7 032 0.9 733 3.893
¿0 • 102, с 1.712 1.468 1.209 1.027 0.5 137
• 104, м2/с 6.528 5.595 4.648 3.917 1.958
т1 • 103, кг 7.121 11.31 20.25 32.97 263.7
а 1.8 2.1 2.55 3 6
Ив 3.221 5.967 12.97 24.85 397.6
Ref 0.2 323 0.3 878 0.8 311 1.667 50.06
•104 1171 632.2 290.8 151.8 9.488
^2•104 4 027 2 536 1416 869.8 108.7
еХ4 • 100 % 0.0833 0.1493 0.3160 0.5799 5.849
е|4•100 % 0.0167 0.0301 0.0617 0.1097 1.058
е^ • 100 % 0.0027 0.0016 0.0096 0.0179 0.343
е?4•100 % 0.0042 0.0075 0.0159 0.0304 0.581
И,е f определялось по формуле
Яе f = Яе (Хтах).
Ошибки в выполнении законов сохранения, возникающие за счет неконсервативности схемы на момент размерного времени Ь = 0.2 с, определялись по формулам
М-2 /-4
^ ] У Лу(ап+2,т — ^п+1,т + сп,т — &п-\,т + е п-2,т)шп,тЩп у/
т=2 п=4
М-2 / -2
У^ Х/(|аП,тШп-2,т| +
т=2 п=2
+ 1 Ьп,,тшп-1,т1 + 1 с3п,тШп,т\ + \(К1тШп+1,т\ + 1 а°п,тип+2,т т
(15)
/ -2 М-4
^ ] ^ ] (ап,т+2 ^п,т+1 + СУ1,т ^п,т-1 + еп,т-2)шп,т^п у/
п=2 т=4
/ -2 М-2
У^ Х/(|ап,тШп,т-2\ +
п=2 т=2
+ \ К,,тшп,т-1\ + \ СУптШп,т \ + Ип,т^п,т+1 \ + \ еп,т^п,т+2 \) Ь
Причем, если вычислялись ошибки в выполнении законов сохранения во всей области, то бралось значение [ = 24 = М, если же они вычислялись для некоторой части жид-
кой зоны, то бралось значение [ = 14. Здесь коэффициенты апт,..., епт, ау
п,т п,т п,т п,т
таковы, что
Л6^п,т
ЬХН 2
п п,т
'(апт^п—2,т ЬптШп-1,т + Спт^п,т ^пт^п+1,т + т^п+2,т)>
Ту ТТ2 (аг1,тШп,т-2 Щг,тШп,т-1 + сп,тШп,т dy,тШ'У,т+1 + еп,т^п,т+2).
Х
3
у
3
у
1
(Множители 1/fa,хп, 1/hут введены для упрощения формул (15).) Для консервативной схемы ер = £Ур = 0.
Из табл. 3 видно, что чем больше фактическое число Рейнольдса, тем больше величины ejg, еУр. Кроме того, эти величины внутри жидкой зоны (0 = 14) значительно меньше аналогичных величин для всей области (0 = 24). Это объясняется тем, что внутри жидкой зоны сетка значительно более мелкая, чем вне ее.
Отметим, что недивергентная форма записи конвективных членов в уравнении для завихренности связана с тем, что фактическое число Рейнольдса Re/ = 50 при г0 = 0.02 м не является малым, так что центральные разности при аппроксимации конвективных членов приведут к пилообразному решению. Второй порядок аппроксимации этих членов выбран для уменьшения влияния схемной вязкости.
На рис. 4 показаны рассчитанные зависимости размерной скорости левого шара от размерного времени для r0 Е {0.006, 0.007, 0.0085,0.01, 0.02} м. Из табл. 3 видно, что безразмерная скорость VL(Kmax) для первых четырех значений г0 меняется слабо. Поэтому размерная скорость УL определяется в основном значением V0:
VL « 0.065V0 = const ■ С2.
Из рис. 4, а видно, что для этих значений г0 скорость левого шара VL с течением времени сначала быстро возрастает, а затем колеблется вокруг приведенного выше значения. Из рис. 4, б следует, что при г0 = 0.02 м скорость левого шара сначала возрастает, а потом начинает убывать. Причем она значительно больше приведенных выше значений, что объясняется большой массой шаров.
На рис. 5, а показаны зависимости расстояния, пройденного левым шаром, от времени для приведенных выше значений г0. Видно, что эти зависимости близки к линейным для первых четырех значений г0. Для г0 = 0.02 м зависимость расстояния, пройденного левым шаром, значительно отличается от линейной.
На рис. 5, б показаны рассчитанные зависимости безразмерного расстояния между центрами шаров от размерного времени. Видно, что чем больше радиусы шаров, тем быстрее они сближаются. Вопрос о значении t, при котором произойдет соударение шаров, остается открытым, так как при г\ ^ 2 шаг г катастрофически уменьшается, не давая величине г\ достигнуть значения 2.01.
Рис. 4. Зависимость скорости движения левого шара от времени при гю = 3 и го = 0.006, 0.007, 0.0085, 0.01 м — кривые 1-4 соответственно (а); г0 = 0.02 м (б)
Рис. 5. Зависимость расстояния, пройденного левым шаром, от времени (а), зависимость безразмерного расстояния между центрами шаров от времени при по = 3 и го = 0.006, 0.007, 0.0085, 0.01, 0.02 м — кривые 1-5 соответственно (б)
-6 -4 -2 0 2 4 6 ^ .6 -4 -2 0 2 4 г -6 -4 -2 0 2 4 г
Рис. 6. Зависимость формы жидкой зоны (темно-серый цвет) вокруг шаров (светло-серый цвет) от времени: £ = 0.257 с (а), £ = 0.9 с (б), £ = 1.8 с (в)
Наконец, на рис. 6, а - в показана форма жидкой зоны при г\0 = 3, г0 = 0.02 м в зависимости от времени ¿. Видно, что она не является симметричной относительно оси г = 0. Эта несимметричность обусловлена значительным вкладом конвективных членов в уравнении для вихря (в 50 раз большим, чем вклад вязких членов), тогда как при их отсутствии задача является симметричной относительно оси г = 0 и скорости шаров одинаковы.
Заключение
Разработан метод расчета нестационарного осесимметричного движения двух твердых шаров в жидкости Бингама под действием силы тяжести. Первым шагом в этом методе является построение подвижной двумерной ортогональной разностной сетки в области течения. Сначала строится дискретный набор так называемых базовых сеток, каждая из которых характеризуется своим расстоянием между центрами шаров Гц, I = 1.15. Далее граничные узлы и некоторые другие характеристики сеток вычисляются с помощью линейной интерполяции соответствующих величин с их ближайших значений для базовых сеток. Затем строится начальное приближение для внутренних узлов с помощью метода трансфинитной интерполяции с границы. Малость расстояния между одноименными узлами текущей и предыдущей итераций по времени обеспечивает
высокую точность начального приближения для построения сетки, так что относительная невязка при использовании всего двух итераций для решении двух уравнений Лапласа лежит в пределах 10-12-10-3. Следующим шагом является решение уравнений Навье — Стокса с вязкостью, зависящей определенным образом от второго инварианта тензора скоростей деформаций в соответствии с регуляризованной моделью жидкости Бингама. Произведены две серии расчетов: первая — для фиксированного радиуса шаров Го = 0.005 м в зависимости от начального расстояния между центрами шаров г 10 Е [2.5, 9] и вторая — для фиксированного безразмерного расстояния между центрами шаров г 10 = 3 в зависимости от радиуса шаров г0 Е [0.006, 0.02] м. Для первой серии расчетов вычислены моменты времени соударения шаров и показано, что чем больше начальное расстояние между центрами шаров, тем медленнее они движутся. Для второй серии расчетов показано, что чем больше масса шаров, тем быстрее они сближаются.
Благодарности. Автор выражает искреннюю признательность В.В. Пухначеву, О.М. Лаврентьевой и В.В. Кузнецову за полезные обсуждения рассматриваемой задачи.
Список литературы / References
[1] Басмат А.С. Нестационарное движение твердых тел цилиндрической и сферической формы в сжимаемой среде: Дис. ... канд. физ.-мат. наук. Киев, 1984. 114 с.
Basmat, A.S. Non-stationary movement of solid bodies of cylindrical and spherical shapes in the compressible medium: Dis. ... cand. of phys.-math. sciences. Kiev, 1984. 114 p. (In Russ.)
[2] Коновалова Н.И. Динамика частиц в вязкой жидкости в быстропеременных полях: Дис. ... канд. физ.-мат. наук. Саранск, 2009. 140 с.
Konovalova, N.I. Dynamics of particles in viscous liquid in the rapidly varied fields: Dis. ... cand. of phys.-math. sciences. Saransk, 2009. 140 p. (In Russ.)
[3] Соковнин О.М., Загоскина Н.В., Загоскин С.Н. Гидродинамика движения сферических частиц, капель и пузырей в неньютоновской жидкости. Аналитические методы исследования // Теор. основы химической технологии. 2012. Т. 46, № 3. С. 243-257. Sokovnin, O.M., Zagoskina, N.V., Zagoskin, S.N. Hydrodynamics of the motion of spherical particles, droplets, and bubbles in a non-newtonian liquid: Analytical methods of investigation // Theor. Foundations of Chemical Engineering. 2012. Vol. 46, No. 3. P. 199-212.
[4] Соковнин О.М., Загоскина Н.В., Загоскин С.Н. Гидродинамика движения сферических частиц, капель и пузырей в неньютоновской жидкости. Экспериментальные исследования // Теоретические основы химической технологии. 2013. Т. 47, № 4. С. 422-433. Sokovnin, O.M., Zagoskina, N.V., Zagoskin, S.N. Hydrodynamics of the movements of spherical particles, drops and bubbles in non-Newtonian liquids. Experimental researches // Theoretical Foundations of Chemical Engineering. 2013. Vol. 47, No. 4. P. 356-367.
[5] Beris, A.N., Tsamopoulas, J.A., Armstrong, R.C., Brown, R.A. Creeping motion of a sphere through a Bingham plastic // J. of Fluid Mechanics. 1985. Vol. 158. P. 219-244.
[6] Пивоваров Ю.В. Расчет движения осесимметричного твердого тела в жидкости Бингама // Вычисл. технологии. 2012. Т. 17, № 4. С. 63-76.
Pivovarov, Yu.V. On calculation of axisymmetric solid body motion in the Bingham medium // Comput. Technologies. 2012. Vol. 17, No. 4. P. 63-76. (In Russ.)
[7] Merkak, O., Jossic, L., Magnin, A. Spheres and interactions between spheres moving at very low velocities in a yield stress fluid //J. Non-Newtonian Fluid Mech. 2006. Vol. 133. P. 99-108.
[8] Liu, B.T., Muller, S.J., Denn, M.M. Interactions of two rigit spheres translating collinearly in creeping flow in a Bingham material //J. Non-Newtonian Fluid Mech. 2003. Vol. 113. P. 49-67.
[9] Malek, J., Ruzicka, M., Shelukhin, V.V. Herschel — Bulkley fluids: existence and regularity of steady flows // Math. Models Methods Appl. Sci. 2005. Vol. 15, No. 12. P. 1845-1861.
[10] Пивоваров Ю.В. Расчет процесса сближения двух сферических капель, находящихся в среде Бингама // ПМТФ. 2014. Т. 55, № 5. С. 29-44.
Pivovarov, Yu.V. Calculating the approach of two spherical droplets located in a Bingham fluid // J. of Applied Mechanics and Technical Physics. 2014. Vol. 55, No. 5. P. 750-763.
[11] Лисейкин В.Д. Алгебраический метод построения разностных сеток: Учеб. пособие. Новосибирск: НГУ, 2002. 46 с.
Liseykin, V.D. Algebraic method for construction of numerical grids. Novosibirsk: NSU, 2002. 46 p. (In Russ.)
[12] Лисейкин В.Д. Метод алгебраической адаптации // Журн. вычисл. математики и ма-тем. физики. 1998. Т. 38, № 10. С. 1692-1709.
Liseikin, V.D. A method of algebraic adaptation // Computational Mathematics and Mathematical Physics. 1998. Vol. 38, No. 10. P. 1624-1640.
[13] Таблицы физических величин / Под ред. И.К. Кикоина. М.: Атомиздат, 1976. 1008 с. Tables of physical quantities / I.K. Kikoin (Ed.). Moscow: Atomizdat, 1976. 1008 p. (In Russ.)
[14] Кочин Н.Е., Кибель И.А., Розе Н.В. Теоретическая гидромеханика. Ч. II. М.: Гос. изд-во физ.-мат. лит., 1963. 727 с.
Kochin, N.E., Kibel, I.A., Roze, N.V. Theoretical hydromechanics. Pt II. Moscow: Gos. izd-vo fiz.-mat. lit., 1963. 727 p. (In Russ.)
[15] Пивоваров Ю.В. Расчет движения жидкости с переменной вязкостью в области с криволинейной границей // Вычисл. технологии. 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.)
Поступила в 'редакцию 1 ноября 2016 г., с доработки —17 января 2017 г.
Calculation of motion of two solid spheres in a Bingham fluid under the action of gravitational force
PIVOVAROV, YURIY V.
Institute of Hydrodynamics of SB RAS, Novosibirsk, 630090, Russia Corresponding author: Pivovarov, Yuriy V., e-mail: pivov@hydro.nsc.ru
The purpose of investigation is the development of algorithm and carrying out the research of flow regimes in axisymmetric non-stationary problem of the motion of two solid spheres in a Bingham fluid under the action of gravitational force.
In the solution of problem the methods of finite differences and the method of alternating directions were emplyed. Solution of vortex problem uses the pentadiagonal matrix algorithm and for the solution of stream function problem so does the Zverev method. In the construction of orthogonal difference grid the alternate triangular method
© ICT SB RAS, 2017
90
K. B. nuBOBapoB
and the method of transfinite interpolation from the boundary were carried out. The derivation of movement equations was done with the help of the method of regularization with a small parameter. The vortex equation was approximated by the convective terms using the upwind differences of the second-order accuracy. Tom's conditions were used as boundary conditions connecting vortex and stream function on the surfaces of spheres. The Reynolds number varies from 0.03 to 50 in the study.
The results of research include the constructed fields of velocity in the flow domain, normal and tangent stresses on the boundaries of spheres, forces acting on the spheres, velocities of spheres, distance between the centers of spheres and the shape of fluid zone around spheres depending on time. For small radii of the spheres, the time of their impact depending on the initial distance between their centers was calculated. It is shown that, for a fixed initial distance between the centers of the spheres, they approach as quicker as their radii are greater.
Conclusion. The mathematical model and numerical algorithm for solving the problem of motion of two solid spheres in a Bingham fluid are developed, the calculations for this model were carried out and qualitative effects of the process of spheres approach were revealed.
Keywords: regularization, non-Newtonian fluid, conformally map, Navier — Stokes equations, vortex, stream function.
Received 1 November 2016 Received in revised form 17 January 2017
