Применение метода R-функций к расчету течений в канавке подшипника Текст научной статьи по специальности «Математика»

УДК 517.95 : 519.63

А. В. АРТЮХ, Н.В. ГИБКИНА, М.В. СИДОРОВ, Е.Е. АГАПОВ

ПРИМЕНЕНИЕ МЕТОДА R-ФУНКЦИЙ К РАСЧЕТУ ТЕЧЕНИЙ В КАНАВКЕ ПОДШИПНИКА

Рассматривается применение метода R-функций для решения задачи расчета плоского стационарного течения вязкой жидкости в канавке подшипника. Вычислительный эксперимент был проведен для единичного квадрата, полученные результаты сравнены с решениями, полученными другими авторами.

Введение

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

Современная вычислительная гидродинамика занимается такими направлениями, как расчет движений вязкой жидкости, численное исследование течений газа с физико-химическими превращениями, изучение распространения ударных волн в различных средах, решение газодинамических задач при наличии излучения и т. д.

В данной работе рассматривается одна из проблем - численное моделирование конвективного течения вязкой несжимаемой жидкости, которое описывается уравнением Навье-Стокса.

Уравнения Навье-Стокса - это система нелинейных уравнений, и построение точных решений является трудновыполнимой, а чаще всего невыполнимой задачей. Обычно точные решения удается получить для простых областей за счет специальным образом подобранной системы координат. В большинстве же своем задачу для уравнений Навье-Стокса приходится решать численно. Будем рассматривать плоские стационарные задачи.

Уравнение Навье-Стокса обычно решают с помощью конечно-разностных (метод сеток) и приближенных аналитических (методы Галеркина, наименьших квадратов и т.д.) методов. Их недостатком является замена сложной области простой приближенной областью. Не свободен от этого недостатка и метод конечных элементов (МКЭ) - один из наиболее употребляемых и эффективных численных методов в механике жидкости. В МКЭ криволинейные участки границы заменяются линейными и при расчете различных сложных конструкций (течения в соплах, форсунках и пр.) приходится измельчать сетку вблизи границы, строить специальные алгоритмы триангуляции и т.п. При построении систем метода конечных элементов обычно применяют либо формулировку метода Галер-кина, либо какой-то вариационный принцип. При этом приходится решать нелинейную систему.

Свободным от этого недостатка является метод R-функций, разрабатываемый школой акад. НАН Украины Рвачева В.Л. Благодаря использованию функций непрерывных аргументов, обладающих свойствами булевых функций, удается строить пучки функций (структуры решений), точно учитывающие геометрию области и краевые условия задачи.

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

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

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

^ ^ 3 Подход смазки

Рис. 1. Схема течения в канавке подшипника: 1 - смазка, 2 - канавка, 3 - подшипник, 4 - вал

При постановке математической задачи делаются следующие допущения:

1. Размеры канавки гораздо меньше диаметра вала, в связи с чем в рассматриваемой области применима декартова система координат.

2. Скорость течения жидкости вдоль канавки настолько мала, что не влияет на структуру решения в поперечном направлении.

3. Ядро холодного масла, медленно текущего вдоль канавки, можно моделировать, задавая низкую постоянную температуру в соответствующей области ядра течения.

Течение масла в канавке подшипника описывается с помощью краевой задачи для системы дифференциальных уравнений относительно функции тока у ( х, у) и температуры Т(х,у) вида (приближение Буссинеска):

дТ

Л2у = Ре1 (Лу, у) + Ог— + РвЬ(х,у) в О, (1)

дх

ЛТ = Ре1 (Т, у) + Рва (х,у) в О, (2)

= £ (8), 8 едО, (3)

= % (8) , ^ дп

Тдо = 1 (8), 8 е дО. (4)

Здесь Ре - число Рейнольдса; Ог - число Грасгофа; Рв - число Пекле; а(х,у) -объемная плотность тепловых источников; 11 (8) - распределение температуры на границе; * ( ) дТх д£у , , а% _

Ь (х, у) =---- ротор вектора внешних сил (I гу ); —, ££ - некоторые распределе-

ду дх 4 7 а8

ния нормальной и касательной составляющих скорости потока соответственно;

л( ч = дЛу ду дЛу ду

J(Лу, у) = ~дх ду ду дх; п - внешняя нормаль к дО. Предполагаем, что Ь, а е Ь2(О), % е'^2(дО), £, 1 е Ь2 (дО), система координат выбрана так, что сила тяжести сонаправлена с вектором (0, -1) .

дО

2. Выбор и обоснование метода решения

Система дифференциальных уравнений Навье-Стокса, описывающая течение масла в канавке подшипника, относительно функции тока у(х,у) и температуры Т(х,у) (приближение Буссинеска) имеет вид:

А 2у = Ре

(

Эх ч АТ = Ре

Ау

Эу

Эу I Эу

Ау

Эу

Эх

Ог — в а,

Эх

( л.Л

Эх

Т ду

--1 Т ^

=% (8), ду

Эп

= § (8), 8 е 5а,

ч Эу | Эу ^ Эх Эу

эп

Т| за = Ь (8) , 8 еэа . Обозначим через (и0, у0 ) решение следующей задачи:

А 2ио = Ог ^ в а,

дх

-Ау0 = 0 в а,

в а ,

эа= %(8) , ^

Эп

= § (8) , 8 еЭа,

'0 |за

= И (8), 8 еЭа.

В задаче (5)-(8) сделаем замену:

у = и0 + и, Т = ^ + ^

где и , V - новые неизвестные функции. Это приводит к задаче

А 2и = Ре

(

Эх

А(и0+и)

Э (и0+и)

Л

Эу I Эу

А(и0+и)

Э (и0 + и)

Эх

_ Эv +Ог— в а,

Эх

-Аv = -Ре

Эх

(Vо + V)

Э(и + и)

Л

и| эп= 0

Эу | Эу Эи

К + V)

Э(и + и)

Эх

в а ,

Эп

= 0,

V Эа= 0.

При замене учитывались условия (9) - (12). В нашем случае краевые условия имеют вид:

[0.5, Эп \{у = 0},

% (8 ) =

0,

у = 0.

% [ 0, ЭП \{у = 0},

§(8) = 1-1, у = 0.

И (8 ) = 1 на ЭП . Структура решения задачи (9) - (12) имеет вид:

(5)

(6)

(7)

(8)

(9) (10) (11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

ЭП

ЭП

и0 = I-ю( +£ ) + ю2Ф0 Уо = Ь + юТо,

_ оЬ д дю д

где =---1---.

дх дх ду ду

Здесь I = ЕС% (8), £ = ЕСЦ (8), Ь = ЕСЬ (8), тогда

0,5у

I = -

у + х(1 - х)(1 - у)

= -х(1- х)(1- у) . у + х(1 - х)(1 - у) Ь = 1.

В структуре решения (20) обозначим

ф = 1-ю(Б11 +£) , а = ю2Ф0, Р = юТо,

тогда и0 = ф + а, у0 = Ь + р . Подставим эти обозначения в (9) - (12):

Л 2 Л 2 ~ Г дЬ др^

Лф+Ла = Ог [дх] в О,

ЛЬ + Лр = 0 в О, да

а = 0,

|дО дп

= 0,

вдО= 0.

Так как в (26) Ь = 1, то (25) - (28) примет вид:

Л2 а = Ог —-Л2 ф в О,

дх

-Лр = 0 в О, да

а = 0,

дО дп

= 0,

Из (30) и (32) получаем, что Р = 0, тогда задача (29) - (30) примет вид:

Л2а = -Л2ф в О,

да

а = 0, —

дО дп

= 0.

Задачу (33) - (34) решаем методом Ритца для уравнения Аа = I. Пусть

а

= ^ СкХ к ,

тогда система Ритца имеет вид:

В нашем случае она имеет вид:

к

Е(АХр Хк ) =(Хк,:1).

и

к

Х(Л% Хк ) = ( -Л 2ф).

и

(20) (21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

(34)

(35)

(36)

дО

дО

дО

Скалярное произведение представляем в виде интегралов:

Е Л-^ХДк^у = -ЦхкА2фdxdy .

]=1 п п

Согласно третьей формуле Грина для оператора Лапласа, которая имеет вид:

-П(Аи - и!Ь (39)

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

ЕАХ)АХкахау = -ЦАХкАФ^у . (40)

]=1 п п

После решения задачи (33) - (34) находим и0, v0 (которое уже известно) и берем их в качестве начального приближения для итерационного процесса.

Рассмотрим теперь задачу (15) - (18). Ее решение также находим с помощью метода Ритца:

к к

и=Е ск тк, v=Е ск Пк.

к =1 к=1

Для (15) система Ритца будет иметь вид:

ЕЕ (а 2 ^ Тк )=

]=1

Тк,Ре

(

Эх

А(и0+и)

Э(и + и)

Эу

Л

| А(и0 + и)

Э(и0+и) Эх

Ог—

Эх

Л

(41)

Или в виде интегралов:

ЕЦА 2 т ]хкахаус]

]=1 П

Л[ Ре

Эх

А(и + и)

Э(и0 + и) Эу

Э(и0 + и)

А(и0+и)

Эу ^ Эх

Ог

Эv Эх

Ткахау. (42)

Интеграл в левой части (42) можно также по третьей формуле Грина (39) представить в виде:

ц а2т ^кахау = цат ^ткахау

Рассмотрим интеграл справа в (42). Используя формулу

г Эу г Эи г

I и—-¿хау = -1 ¿хау + I uvco8(n,x)dxdy

П Эх П Эх 1

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

(43)

(44)

111 Ре

-А(и + и)

Э(и0 + и) ЭТк Эу Эх

" А(и + и)

Э(и0 + и) ЭТк Эх Эу

Ог ^

Ог—т

Эх

dxdy

или после преобразований

III РеА(и0 + и)

Г Э(и0 + и) ЭТк Э(и0 + и) ЭТк

Эх Эу Эу Эх

Ог—т

Эх

dxdy

(45)

Для (16) система Ритца имеет вид:

-Ре

Эх

(vо + ^

к

Е(-АПр Пк )

]=1

Э(и0 + и) Эу

Э(и + и)

ЭуI ( 0 ) Эх

(46) 29

п

п

п

Или в интегральной форме:

-Ж*

Sx

(vo +v)

Z| -ЛАПA^dy

j=1 V n

d(Uo + u) Sy

cj =

,_d_f (v +v) d(Uo + u)

Sy V 0 dx

Hkdxdy.

(47)

Интеграл в левой части (47) можно преобразовать по первой формуле Грина для оператора Лапласа:

-[ vAudxdy = jf —U— + —— dxdy - j vdudS

n nVdx dx Sy SyJ S dn

Тогда этот интеграл будет равен:

-jjAnJnkdxdy = jj|

'^nj dnk ^nj on.

Л

dxdy

J

О О[дх дх ду ду

Интеграл в правой части (47) также преобразуем по формуле (44):

(48)

(49)

-jj[pe

-(vo + v)

д(uo + u) дПк

Sy Sx

+ (vo + v)

d(uo + u) дПк

Sx Sy

dxdy .

Или после преобразований:

■jl|Pe

(vo + v)

d(uo + u) дПк d(uo + u) дПк

Sx

Sx

dxdy

(5o)

ду ду

Неопределенные компоненты в (20), (21), (23), (24) аппроксимируем кубическими сплайнами Шенберга.

3. Результаты численного эксперимента

Вычислительный эксперимент был проведен для прямоугольной области. Везде предполагалось, что массовые силы Г потенциальны, т.е. гОГ = 0, а значит, Ь(х,у) = 0, а (х,у) 0 . В качестве базисных функций выбирались степенные полиномы, тригонометрические полиномы, полиномы Лежандра, кубические сплайны Шенберга. Шаг сетки сплайнов выбран Ьх = Ьу = 0,1. При вычислении интегралов в системах Ритца использовалась формула Гаусса с 16 узлами по каждой переменной на каждом частичном квадрате 0,1 х 0,1. Предполагалось, что жидкость ограничена твердыми неподвижными стенками. Начальное приближение было построено для чисел Рейнольдса Ре = 1, Пекле Ре = 1 и Грасгофа Ог = 50 .

Вычислительный эксперимент был проведен при а = Ь = 1. Результаты в виде линий уровня температуры (изотерм) и линий уровня функции тока (линий тока) приведены на рис. 2 (Ре = 1, Ре = 1 и вг = 50 ).

1-"

0.6

0.4

0.2

0 0.2 0.4 0.6 0.S 1 0 0.2 С. 4 0.0 С.(

а б

Рис. 2. Изотермы T = const (а) и линии тока у = const (б)

n

n

Полученные результаты хорошо согласуются с результатами физических экспериментов [1, 2] и результатами, полученными методом фиктивных областей [9, 10]. Расхождения составили около 5%.

Выводы

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

Литература: 1. ЛандауЛ.Ф., ЛифшицЕ.М. Теоретическая физика. В 10 т. Т. VI. Гидродинамика. М.: Физматлит, 2003. 736 с. 2. ЛойцянскийЛ.Г. Механика жидкости и газа. М.: Дрофа, 2003. 840 с. 3. Кочин Н.Е., КибельИ.А., РозеН.В. Теоретическая гидромеханика. Ч. 2. М.: ГИФМЛ, 1963. 728 с. 4. ЛамбГ. Гидродинамика. М.: РХД, 2003. Т. 1. 452 с.; Т. 2, 452 с. 5. ШкадовВ.Я., Запрянов З.Д. Течения вязкой жидкости. М.: Изд-во Моск. у-та, 1984. 200 с. 6. Свешников А.Г., Боголюбов А.Н., Кравцов В.В. Лекции по математической физике. М.: Изд-во МГУ; Наука, 2004. 416 с. 7. Михлин С.Г. Вариационные методы в математической физике. М.: Наука, 1970. 512 с. 8. Рвачев В.Л. Теория Я-функций и некоторые ее приложения. К.: Наук. думка, 1982. 552 с. 9. Вабищевич П.Н., Вабищевич Т.Н. Численное решение стационарных задач вязкой несжимаемой жидкости на основе метода фиктивных областей // Вычислительная математика и математическое обеспечение ЭВМ. М.: Изд -во Моск. ун-та, 1985. С. 255-262. 10. Вабищевич П.Н. Метод фиктивных областей в математической физике. М.: Изд -во Моск. ун-та, 1991. 156 с.

Поступила в редколлегию 22.06.2009 Артюх Антон Владимирович, студент группы ПМм-08-1 факультета ПММ ХНУРЭ. Научные интересы: математическое моделирование, математическая физика, численные методы. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. 702-14-36.

Гибкина Надежда Валентиновна, канд. техн. наук, доц. каф. прикладной математики ХНУРЭ. Научные интересы: экономический риск, актуарная математика, математическая физика, оптимальное управление динамическими объектами. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. 702-14-36.

Сидоров Максим Викторович, канд. физ.-мат. наук, ст. преп. каф. прикладной математики ХНУРЭ. Научные интересы: математическое моделирование, математическая физика, теория Я-функций и ее приложения. Увлечения и хобби: история культуры. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. 702-14-36.

Агапов Евгений Евгеньевич, студент группы ИНФ-05-2 факультета ПММ ХНУРЭ. Научные интересы: численные методы, компьютерное моделирование. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. 702-14-36.