CC BY
48
Вычислительные технологии
Том 18, № 2, 2013
Численный расчёт свободного движения малого
•• U U
объёма вязкой несжимаемой жидкости между вращающимися цилиндрами
А. В. Плничкин1, Л. Г. Варепо2 1 Омский филиал Института математики им. С.Л. Соболева СО РАН, Россия 2Московский государственный университет печати, Россия e-mail: [email protected], [email protected]
Проведено моделирование течения, и исследована начальная картина протекания малого объёма вязкой несжимаемой жидкости, имеющей свободные границы, между вращающимися цилиндрами при числах Рейнольдса от 1 до 100 на двумерной сетке с помощью конечно-разностных методов с равномерным шагом. Получены численные решения движения свободных границ жидкости, характерные для различных значений вязкости, путём расчётов перемещения граничных узлов по узловым линиям фиксированной сетки.
Ключевые слова: моделирование, несжимаемая вязкая жидкость, число Рейнольдса, конечно-разностная сетка, уравнения Навье — Стокса.
1. Постановка задачи
Рассмотривается моделирование процесса переноса жидкости на основе решения уравнений Навье — Стокса для вязкой несжимаемой жидкости. Поскольку анализируется движение жидкости с начальным нанесённым её слоем на верхний цилиндр, целесообразно рассмотреть область решения в сопутствующей системе координат ОХУ, связанной с верхним цилиндром с центром в точке О, расположенной на поверхности цилиндра в середине нанесённого участка жидкости (рис. 1).
Ограничимся постановкой плоской задачи в полярной системе координат г, в, связанной с центром О1 цилиндра 1 и вращающейся с постоянной угловой скоростью ш в положительном направлении ((р = в + ш£). При этом отсчёт углов производится от отрицательного направления оси ОУ. Уравнения Навье — Стокса для вязкой несжимаемой жидкости в неподвижной полярной системе координат г, для вектора скорости (У, Ур) имеют вид ([1], с. 363)
ЗУГ
"at
+ vr
dVr
dr
r V,R dVr
r
d'R
(V,)2
1 dP
p dr
—+ И v2Vr -
Vr 2R 3Vin
r2 r2 d'R
dt
, + V dVk + VvRdVV^ + VVv
+ Vr Г\ + Г\ -F-, +
dr
r
d'R
r
1 R dP p r d'R
+ И V2V,-^ +
V, 2R dVr
r2 r2 d'R
dVr Vr R dV, dr r r d'R
r
где
V2 V = -—+
1 дУ д2 V Я2 д2 V
+
г дг дг2 г2 (д^Я)2' V — кинематическая вязкость, р — плотность жидкости, Р — давление, Я — радиус цилиндра 1.
При переходе к сопутствующей полярной системе координат г, в будем использовать уравнения для относительных компонент скоростей Уг,Уе. Перевод компонент вектора ускорения жидкости в эту систему координат с угловым ускорением г и угловой скоростью ш без изменения начала координат выражается уравнением
дУг дУ\
дг' дг
дЦ- д~и) + ( 0,гг ) + ( —2шЦ, 2шЦ ) + ( —ш2 г, 0
где предпоследнее слагаемое в правой части представляет кориолисово, а последнее — центростремительное ускорение. После перехода к новым компонентам скорости с преобразованием ускорений и при учёте углового ускорения г рассматриваемые уравнения Навье — Стокса примут следующий вид:
дит
дг
+ Ц
дит
т ЦЯ дит (Цв + иг)2 1 дР ( 2 ит 2Ядив
- — +--— - —---+ VI V и - — —7Г — —й--
дг г двЯ г р дг \ г2 г2 двЯ
див
дг
+ и-
див ив Я див и- ие
дг
г двЯ г
+ 2итш+гг — —
1 Я дР
ди- и-
р г двЯ Ядив
Ы V2 ив
ц 2ЯдиЛ Т2+Т2 двЯ)
дг + г + г двЯ
— 0,
(1)
, (2) (3)
где
2и —
1 ди д2 и Я2 д2 и
+
+
г дг дг2 г2 (двЯ)
2
Цилиндры 1, 2 имеют радиус R и вращаются в противоположные стороны с угловой скоростью ш без ускорения (е = 0), при этом наименьшее расстояние между поверхностями цилиндров равно 8. Рассмотрим прикасающуюся к цилиндру 1 область жидкости П толщиной по радиусу 8s , ограниченную частью окружности длиной 8l (8 < 8s). В начальный момент времени скорость жидкости в сопутствующей системе координат r, в равна нулю, на свободных поверхностях, где жидкость не соприкасается с границами цилиндров, P = Ратм (Ратм = 105 Н/м2). Для перехода к безразмерным переменным введём характерные длину L = 8s и скорость V0 = шR, число Рейнольдса Re = V0L/v, число Эйлера Eu = PaTM/pV0 = P0. Эти величины будут использованы в качестве параметров в тестовых расчётах для сравнительных характеристик.
Координаты точек окружности цилиндра 2 с центром в точке O2, расположенной на расстоянии (2R + 8) от O1, представим через (r, ф-\), где ф1 — угол между направлением из O1 на рассматриваемую точку цилиндра 2 и отрезком, соединяющим O1 и O2. Эти координаты имеют вид
r = Ryjsin2 Ф2 + (2 + 8/R - cos)2 = R^5 + 48/R + 82/R2 - 2(2 + 8/R) cos^, (4)
= arcsin(sin^2/V5 + 48/R + 82/R2 - 2(2 + 8/R) cosФ2), (5)
где ф2 — угол соответствующей дуги на цилиндре 2.
В сопутствующей системе координат движение центра цилиндра 2 вокруг цилиндра 1 происходит с угловой скоростью ш, а вращение его вокруг своего центра — с угловой скоростью 2ш по часовой стрелке. При соприкосновении цилиндра 2 с областью жидкости учитывается скорость движения точек окружности цилиндра 2, имеющая следующие компоненты в системе координат r, в:
V = Vx cOs('m1) + Vy sin('m1) , Vr = Vx sin('m1) - Vy cOs('m1), (6)
для которых xm = Xc + (2R + 8) sin('1) - r2 sin('2), Ут = Ус - (2R + 8) cOs('1) + r2 cOs('2), Vx = -ш(2R + 8)cos('1) - 2шг2 cos('2), Vy = -ш(2R + 8) sin('1) + 2шг2 sin('2), '1 = '0 - шí, '2 = 'mo - 2ш£; (xm, ym), (xc, yc) — координаты рассматриваемой точки и центра цилиндра 1 в декартовой системе координат OXY, связанной с цилиндром 1 (см. рис. 1), где '1 — угол поворота центра цилиндра 2 O2 относительно центра цилиндра 1 O1 (угол в для O2), '2 — угол между направлением из O2 на рассматриваемую точку цилиндра 2 и вертикальной линией, '0 — начальное значение в для центра цилиндра 2, 'm0 — начальное значение '2.
С изменением со временем положений точек xm,ym по углам '1, '2 их координаты в полярной системе координат r, в определяются следующим образом. Сначала по (4) для каждого значения угла ф2 = '2 - '1 находится r, затем по (5) вычисляется ф-\_. Окончательно имеем в = 'm1 = '1 -ф]_. При этом радиусы на поверхностях цилиндров 1 и 2 равны r1 = r2 = R.
Для свободной границы жидкости П, представляемой в виде некоторой функции f (t, r, в), кинематическое условие, обеспечивающее непроницаемость границы, будет иметь вид
ft + fr Ur + fe Ue = 0. (7)
Сила поверхностного натяжения жидкости на свободных границах для цилиндрической поверхности определяется по отношению Cn/rcr, где Сп — коэффициент поверхностного натяжения жидкости, rcr — радиус кривизны линии f (t, r, в).
При слабом взаимодействии свободной границы жидкости с окружающим газом, имеющим давление Ратм и малые плотность и вязкость, условие непрерывности тензора напряжений на границе двух сред можно записать в виде формулы Лапласа
Р = Ратм — Сп/т"^гГсг ' п (8)
где п — внешняя нормаль к границе.
Для области жидкости П в сопутствующей системе координат на начальный момент времени Ь = 0 координаты т, в имеют пределы т Е [Я, Я+8в], в Е [—8ь/(2Я), 8ь/(2Я)]. Во всей области иг(0, т, в) = 0, Ид(0, т, в) = 0. В последующие моменты времени цилиндр 2 станет соприкасаться с этой областью начиная с точки (Я + 8s,8ь/(2Я)), для которой из условий прилипания и непротекания компоненты скорости будут определяться по формулам (6) иг = У, Ид = Уд. Это же будет выполняться и для других граничных точек жидкости, приходящих в соприкосновение с подвижной границей цилиндра 2. Условие для градиента давления на этой границе для т, в, соответствующих значениям (4) и (5) при ф\ = — в, получается из уравнений (1), (2):
(1 дР . лч Я дР . ,Л ((Ид + шт)2 Г2тт иг 2Я дИд \ (рдРР(Ь,т,в),7тдш(Ь,т,в)) = \;—- + 'Г*— й — йШ) ,
Иг Ид птт , Ид , 2ЯдИг
- 2ИГш + V Ч2ид- -д +
т2 т2 двЯ
На границе с цилиндром 1 из тех же условий для Ь из [0,Т] граничные условия в сопутствующей с этим цилиндром системе координат примут вид
1дР д 2П
иг(г,я,в) = 0,Ид(ь,я,в) = 0, -—(ь,я,в) = ш2я + и—т(ь,я,в).
р дт дт2
При этом координата в будет изменяться в пределах соприкасания рассматриваемой области жидкости с цилиндром 1.
На свободной границе необходимо поставить условия для компонент скорости, для каждого момента времени определяемые из следующих условий согласования тензоров напряжений на свободной границе (отсутствие взаимодействия с внешней средой): дит/дп = 0, дип/дп = 0, где ит — тангенсальный вектор скорости на границе жидкости, ип — вектор скорости по нормали к границе. При этом давление на свободных границах для всех Ь из [0,Т] определяется из упрощённого динамического условия (8).
При использовании уравнения (3) для определении поля давления во всей расчётной области течения жидкости П его обычно дополняют эволюционным членом дР/дЬ в виде
е,дР + V ■ и = 0, (9)
где
и = ^ )• V- и = ^ + И + ^
дт т т двЯ
ер > 0 — параметр, оптимально выбираемый при численных расчётах для сходимости решения и приближения (9) к уравнению (3).
2. Конечно-разностные методы решения
Для численного решения уравнений (1), (2), (9) вводится новая полярная система координат х,у, преобразованная из полярной системы координат г, в и соответствующая (Яв, Я — г), в которой граница цилиндра 1 проходит по оси Ох, а у цилиндра 2 центр перемещается при у < 0 от Яр0 справа налево со скоростью Яш. В этой системе координат расчётную область W представим в форме прямоугольника с регулярной сеткой и равномерными шагами Нх, Ну (Мх, Ыу — число узлов по координатам х, у). На фиксированной сетке применяются конечно-разностные методы с вводом подвижных граничных узлов для границы цилиндра 2 и свободной границы жидкости, которая в начальный момент находится на цилиндре 1 без относительного движения. В новой системе координат обозначим компоненты скорости (Ц, — Ц) через (и, у) = V.
В принятых в системе координат х,у обозначениях векторов скорости (что удобно для привязки к размерной длине области жидкости при графическом отображении в прямоугольной расчётной области) уравнения (1)-(3) после перестановки (1) и (2) и при отсутствии углового ускорения е примут следующий вид:
ди ди иЯ ди уи ^ дЬ ду Я — у дх Я — у
1 Я дР / 2 и 2Я ду \ , Л
+ А У2и-—---- , (10)
р Я — у дх \ (Я — у)2 (Я — уг2 дх) '
ду ду иЯ ду (и + ш(Я — у))2 1 дР / 2 у 2Я ди
т +%+ЯГ-уь + я-у = — ррду П у2у — (Я—у)* Э1-.
ду у = 0, (12)
где
ду Я — у Я — у дх
„2тт 1 ди д2 и Я2 д 2и V2U = —--+ +
Я — уду ду2 (Я — у)2 (дх)2'
Для расчёта в уравнениях (10) и (11) конвективно-диффузионных членов использовалась схема стабилизирующей поправки [2] с итерационным шагом по времени т
Vn+1/2 — V
Л I \ / П+1 / 2 Ж / П \ | \ \ / П |
т
Лl(vn+1/2 — Vn) + ЛVn + Л0Vn — Грп+1/р + ¥п, (13)
^^п+1 _ vn+1/2
- = л2
Л2(Vn+1 — Vn), (14)
где
л РЯР ^ + Л_1 2КЪ \2 А1А-1 л А2 + А_2 , А2А-2
Л1 = иЯ/(Я — у)^Х" + иЯ/(Я — у) -АТ, Л2 = у + '
ЛoVn = А —2Я/(Я — у)2 -+Л-■ уп — 1/(Я — у) -2hЛ_2 ип—
у
— 1/(Я — у)2пп, 2Я/ (Я — у)2 ^ ^ — 1/(Я — у)А2 +ьА_2 уп — 1/(Я — у)2уА +
ху
+ ^1/(Я — у)ьпип, —1/(Я — у)ипип^
Еп = (^2упш, —2ипш — ш2(Я — у)
А\, Д-1, Д2, А-2 — операторы сдвига функции на шаг сетки вверх или вниз по осям координат х, у.
Оператор Л в (13) является суммой Л1+Л2. Порядок аппроксимации в операторах Л и Л0 не больше двух (т.е. 0(кХ,Щ)), а порядок аппроксимации градиента давления (Грп+1 в (13)) и дивергенции скорости в (9) в случае применения трехточечных шаблонов в каждом пространственном направлении может быть повышен до четырёх. Для градиента давления используется следующее математическое представление в окрестности узла (хг}уу):
ТРп+1| ^ (Я/(Я у)(Д1 + Д-1 Рп+1 К +1 N Д2 + Д-2 рп+1 Ч гп+1 \ (15
ГР кз - ^Я/(Я — у) ^ 2ьх — ~6рхххм) , 2ку — ~6РУУУ^) . (15
Для аппроксимации третьих производных по x и y от давления второго порядка на компактном шаблоне (с использованием трёх узлов в каждом направлении) вначале произведём замены этих производных в виде
Pxxx = ((R - y)2/R2) (pGx - Pyyx + R/(R - y)pyx), (16)
Pyyy = pGy - 2R2/(R - y)3pxx - R2/(R - y)2Pxxy + 1/(R - y)2Py + 1/(R - y)Pyy, (17)
где G следует из подстановки dUr/dt и OUg/dt из уравнений (1) и (2) в (3) и равняется всем слагаемым без функции давления после преобразования координат и компонент скорости ((9R,R - r) в (x,y) и (U,-Ur) в (u,v)):
G(x, y, u, v) = -v2y + 1/(R - y)(vvy - 2uuy - 2ruyvx) +
+R/(R - y)2(uvx + uxv - Ru2x) + 2ш2 - 2шщ + 2w/(R - y)(u + Rvx). (18)
Аналогично строится аппроксимация дивергенции скорости четвёртого порядка на компактном шаблоне
V ■ V"+1|„ ^ R/(R - y)( Al +¡A-1 - f K+Lj
-1/(R - vKf + v"+ - f v"»., (19)
с использованием замен для третьих производных и с последующей аппроксимацией первых и вторых производных по каждому пространственному направлению центральными разностями второго порядка
иххх = —(Я — у)/ЯУухх + 1/ЯУхх, (20)
Уууу = 2/(Я — у)3(у — Яих) + 2/(Я — у)2(Уу — Яиху) + 1/(Я — у)(Ууу — Яихуу). (21)
Аппроксимация производных в (15) и (19) производится в узлах с индексами г,], где г = 1,..., N — 1, ] = 1,..., N — 1. Для расчёта давления уравнение (9) представляется в конечно-разностном виде
ерр-Р-+ V • Уп+1,к+1 = 0, (22)
т
где дивергенция скорости заменяется по (19) при использовании (20), (21) с заменой частных производных первого и второго порядка на конечно-разностные аналоги со вторым порядком аппроксимации на трехточечных шаблонах. При каждом расчёте давления необходим дополнительный расчёт данного параметра на твёрдых границах по уравнениям движения (10) и (11). Для этого вблизи твёрдых границ в дополнительных узлах на расстоянии в полшага от стенок используется аналог условия Тома для производных от компонент скорости и для их значений в виде (на примере граничного узла (0,]))
д'2 « \ = «Ц — „ = (31'0„,- + б1'ц- — У2,з)
(сдх)2) у ' 1/2' 8 .
На каждой п + 1-й итерации по времени по (22) определяется давление с помощью отдельного итерационного процесса при к = 0,1, ...,Ж(е), где £ — малая величина, не превышающая по норме приращения давления на Ж-й итерации. При этом на каждой итерации по к проводится перерасчёт вектора скорости уп+1,к при изменённом давлении рп+1,к. Данный алгоритм был применён в работе [3] для тестового расчёта вязкой несжимаемой жидкости. В такой постановке моделирование течения вязкой несжимаемой жидкости производится до определённого значения времени £, пока цилиндр 2 находится в зоне контакта с жидкостью, первоначально расположенной на цилиндре 1 и при начальном времени вступающей в контакт с цилиндром 2.
3. Расчёт движения границы жидкости
Для расчёта перемещений условной границы жидкости, не проходящей через узлы сетки, вводятся дополнительные узлы на линиях между внутренними и внешними узлами, как это показано на рис. 2. Их смещение за шаг по времени т с учётом с каким-либо приближением кривизны границы можно определить по компонентам скорости («^ )
Рис. 2. Сеточный шаблон около подвижной границы
в прилегающем к границе расчётном (г, ])-м узле. Для этой цели можно пронумеровать окрестные узлы от 0 (центральный (г, ])-й узел) до 8 и произвести построение параболической интерполяции для кривой границы по трём дополнительным узлам, находящимся на пересечениях с узловыми линиями. В таком случае разрешение малых структур для областей, занимаемых жидкостью, будет ограничиваться шагами по сетке Нх и Ну. Как и центральный узел, указанные три узла выбираются на данной узловой и двух соседних линиях.
Из рис. 2 видно, что для узловой линии такие узлы по оси х имеют номера 36, 02 и 04, по оси у — 02, 04, 17. Пусть при рассмотрении дополнительного узла по оси х будут заданы координаты этих узлов (х1,у1), (х2,у2), (х3,уз). Тогда интерполяцию можно записать в следующем виде:
х = Х2 + (у — У2) (Х3—Х1) + (*—^ (хз—Х2 — . (23)
V Уз — у1 / уз — у1 V уз — у2 у2 — у1 )
Аналогично интерполяция строится и для дополнительного узла по оси у. В новых координатах такие интерполянты вместо f (¿, г, в) для кинематического условия (7) можно рассматривать как отдельные функции свободной границы f1(t,x,y) и ^^,х,у), причём одна из них представляет перемещение свободной границы по оси х, другая — по оси у.
После перемещения границы в виде интерполянты (23) за время т координаты рассматриваемых узлов по х и у изменятся на величины тщ^ и ту^^ и примут значения (х'1,у'1), (х2,у2), (х'3,у3). При этом интерполянта для узловой линии по оси х для определения новых координат дополнительного узла (х2,у2) сместится на величину
гр^ _ А^^ \ (^г/1) 2 / т^ _ /у»? гр? _ гр?
5х = тщ,! — туи( + ^^ ( — ) . (24)
. уЗ — у 1 / уЗ — у1 V уЗ — у2 у2 — у 1
Представленный расчёт движения границ жидкости по дополнительным узлам, находящимся на узловых линиях, позволяет учесть кривизну линии границы и точнее вычислять смещение последней при её произвольных ориентациях и произвольных направлениях вектора скорости.
4. Результаты расчётов
Для рассмотренной постановки задачи приведём результаты расчётов по конечно-разностной схеме (13), (14), (22) на равномерной сетке при Мх = Му = 80 и V в пределах от 0.2 • 10-5 до 0.2 • 10-3. Для определения числа Рейнольдса используем такие параметры как начальная толщина слоя жидкости и начальная радиальная скорость схождения цилиндров в области жидкости. При выбранных значениях этих параметров течение будет соответствовать числам И,е от 1 до 100, причём скорость между цилиндрами в малой области их контакта с жидкостью по мере вращения цилиндров может изменяться в несколько раз.
Для итерационного шага т значения выбирались из условий устойчивости расчётной схемы в пределах от 0.5 • 10-7 до 0.5 • 10-6 при следующих параметрах задачи: начальные размеры области жидкости 8ь = 0.008 и 8в = 0.004, коэффициент поверхностного натяжения Сп = 0.03, г1 = г2 = 0.05. Размерности всех величин здесь соответствуют
д
Рис. 3. Границы области течения жидкости при г = 0.1 ■ 10-2, Яе = 1 (а), г = 0.1 ■ 10 2, Яе = 10 (б), г = 0.2 ■ 10-2, Яе = 10 (в), г = 0.1 ■ 10-2, Яе = 100 (г) и г = 0.2 ■ 10-2, Яе = 100 (д); ш = 100
размерностям физических величин в системе СИ. Расчётная сетка рассматривалась для области, включающей начальную область жидкости П, с координатами х от —0.008 до 0.008 и у от —0.0049 до 0.0001.
На рис. 3 показаны численные решения с мгновенными линиями тока, полученными при задании на одной границе ф = 0 из поля скоростей путём интегрирования в расчётной области с точностью 0(Н2х + Ну) следующих соотношений:
я—удф = щ дф = (25)
Я ду ' дх '
Визуализация течения жидкости между двумя вращающимися цилиндрами, представленная на рис. 3, показывает возможные перемещения свободных границ жидкости для разных чисел И,е на моменты времени t = 0.1 • 10-2 и t = 0.2 • 10-2. При этом движение свободных границ при числах И,е = 1 и 10 по сравнению с расчётами для И,е = 100 соответствует течению с большей вязкостью и существенным изменением всей свободной границы на начальном промежутке времени движения цилиндров.
Таким образом, показаны возможности применения рассмотренного метода для моделирования течений несжимаемой жидкости с движущимися и свободными границами на заданной регулярной сетке.
Список литературы
[1] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1978. 736 с.
[2] Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967. 196 с.
[3] Плничкин А.В. Ускорение сходимости в расчётах стационарных течений жидкости при больших числах Рейнольдса // Вычисл. технологии. 2008. Т. 13. Спец. выпуск № 3. С. 38-44.
Поступила в 'редакцию 7 декабря 2012 г., с доработки — 25 февраля 2013 г.
