Численное решение уравнений движения методом конечных элементов в условиях ламинарного течения вязкой жидкости в конвергентном канале Текст научной статьи по специальности «Физика»

ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ ДВИЖЕНИЯ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ В УСЛОВИЯХ ЛАМИНАРНОГО ТЕЧЕНИЯ ВЯЗКОЙ ЖИДКОСТИ В КОНВЕРГЕНТНОМ КАНАЛЕ

А.Г. БАГОУТДИНОВА, Я.Д. ЗОЛОТОНОСОВ

Казанский государственный энергетический университет

В работе описан алгоритм решения трехмерной задачи течения вязкой несжимаемой жидкости в криволинейном конвергентном канале методом конечных элементов.

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

В связи с этим возникает необходимость создания математической модели и алгоритма численной реализации задачи гидродинамики ламинарного течения вязкой жидкости в таких каналах.

Рассмотрим течение жидкости в криволинейном конвергентном канале, вращающемся вокруг своей оси, перпендикулярной направлению движения

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

поверхность вращения, заданную уравнением г = к(г), которая снабжена радиальными лопатками. На выходе из канала расположена кольцевая насадка с ромбовидными отверстиями (рис.1, 2).

Рис. 1. Конвергентный канал

© А.Г. Багоутдинова, Я.Д. Золотоносов Проблемы энергетики, 2007, № 3-4

Рис. 2. Сечение канала

Обозначим через у = в — а угол, образованный соседними лопатками, п -число отверстий, к - номер отверстия. Тогда “ширина” й1 одного отверстия будет

уЯ / ч

равна ----, а “высота” й2 будет равна к(Я). Область отверстий может быть

п

записана следующим образом:

2h(R)(пфі у(k 11)) у (k 11) у (k ~ 1/2)

-------------------, при--------<ф<-------------

D: Id <

у n n

2 h(R )(ky і пф) у (k і 1/2 ) yk

-----------------, при-------------<ф<------

у n n

, k = 1,2...,n ,

а область жесткой стенки вокруг этих отверстий так:

’ 2h(R)(пф~ у (k 11)) у (k 11) у (k ~ 1/2)

-------------------, при---------<ф<-------------

D:\z\>

у n

2h(R)(ky і пф) у (k і 1/2)

---------------, при------------

у n

<ф<

yk

, k = 1,2..., n

Запишем уравнения движения и неразрывности для течения вязкой жидкости в конвергентном канале в цилиндрической системе координат г, ф, г с учетом центробежных и кориолисовых сил [1, 3]:

r ф ^' r —+—-----

1 d p

=--------+ v

p d r

ґ

AVr і

dVr T ' i 2

Z dz r

2 dV9 Vr'

r 2 дф " r 2.

+ 2ю(r sinф + Vф cos ф)+ ю 2r cos ф ; (1)

SV9 VpSVp дУф VrV9

Vr——+———+Vz——+-----—

dr r дф dz r

n

n

r

1 дp =----------------------+ V

рг дф

ЛУф+" 2-

2дVr Vф "1

г дф г дУг УфдУг дVz 1 д p

— 2ю(Уг сое ф — Уф «іпф++ ю2г «іпф ;

-+У,

=--------+ V АУг;

дг гдф дг р дг

(2)

(3)

дУг дУг Уг 1 дУф

—+—+—+------- =0;

дг дг гг дф

(4)

д2 1 д 1 д2 д2

- + -

где Д = —- +---------------------------------------+ - 2 2

дг г дг г дф дг

- оператор Лапласа.

В качестве условий однозначности для системы уравнений (1)-(4) зададим профили скорости, давления во входном сечении канала и граничные условия на стенках:

на входе (г = го):

Уг = ио, Уф = °> Уг = °; р = ро;

на стенках канала:

(ф = а) : Уг = 0, Уф = ю г, Уг = 0; (ф = в) : Уг = 0, Уф = ю г, Уг = 0;

(г = Н) : Уг = 0, Уф = юг, Уг = 0;

(5)

др дУг

в плоскости симметрии (г = 0) : Уг = 0, — = 0, ----------------

дг дг

= 0,

дУф

дг

= 0;

на выходе (г = Я)

через отверстие: Уг = иср, Уф = юЯ, Уг = 0;

на жесткой стенке: Уг = 0, Уф = юЯ, Уг = 0.

Здесь иср = О-, где 2 - расход жидкости в канале, £ - площадь отверстия. пБ

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

дУг дУг Уг 1 дУф р

-------+-------+ — +-----------— =------,

дг гг дф к

дг

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

р = — к

ґ дУг дУг Уг 1 дУф 4

------+-----------+ — +-------—

^ дг дг гг дф

Подставив полученное выражение для давления в уравнения движения (1)-(3), получим:

ЗУг Уф ЗУг дУг Уф к

Уг —— +^- —+Уг —— —5- = -

дг г дф дг г р

З 2УГ д 2Уг 1 дУг Уг 1 д 2Уф 1 дУф

—-—\-----------------------\-—+------- ------

дг дгдг г дг г г ЗфЗг г дф

+

+ V

ґ 2 дУф Уг Л

ДУг---------------Ї---

г 2 дф г

\

2

дУф Уф дУф дУф УгУф к

Уг——+——— + Уг——+------— = —

дг г Зф дг г р г

+ 2юУ «іп ф + Уф со« ф)+ ю2г со«ф;

ґд 2Уг д 2У7 1 ЗУг 1 д 2У 1

ф

Зг дф дг дф г дф г дф2

(8)

+

+ V

АУ_ +

2 дУг Уф

ф 2

г дф г "

дУг УфдУг дУг к

У^—- + ф г + У^—1 = -Зг гдф дг р

— 2ю(уг со« ф — Уф «іпф)+ ю2г «іпф; ґд 2Уг д 2Уг 1 дУг 1 з 2 У 1

г г и* г

— + — + +--

ф

дг Зг Зг2 г дг г дфдг

+ vAyz.

(9)

(10)

Решение системы (7)-(10) будем искать в виде

Уг = и0 / (г.ф ,г ), Уф = югО (г,ф ,г ),

Уг = и0 Н (г,ф ,г), Р — Р0 = и0 р Р (г,ф ,г +

(11)

_ г _ г

Здесь г =—, г = — - безразмерные переменные; и о - начальная Го го

скорость; ро - давление на свободной границе; /,С,Н,Р - соответственно безразмерные компоненты радиальной, тангенциальной, осевой скоростей и параметра давления.

Подставляя вид решения (11) в систему (7)-(10), получим:

д/ д/ д/ 2 2 /— + Ю— + Н— — 12О 2 г = к

Зг Зф зг

Гд2 / З2 Н З2 в / 1 З/Л

—— + —— +1 —— + —

дг Згдг ЗфЗ г г г Зг

+

1 Г 2 Г дО /

ш

Ие [ г дф г2

.дО дв дО

+ ГОР— + рн

дг дф дг

■ + 2Г/зіп ф +Г (2в + 1)г еозф ;

ГР

( д2/ д2Н д20 1 д/Л

_— +---------------+1— + _—

дгдф дгдф дф2 г дф

+

(12)

1 Г де 2 д/]

+ —^ гДв + 2-------+-----------^- 2/еозф +1(20 + Р)зіпф ;

Ие [ дг ГР2 дф ]

(13)

дН дН дН

/----------+ Ю----------+ Н = 1

дг дф дг

Гд2/ д 2Н д 1 д/Л

+ Г

-+-

-+—

дгдг дг2 дфдг г дг

1

+—ДН; Ие

(14)

(д / дН дО /'

Р= -1 —+----------+ Г—+—

дг дг дф г ,

(15)

с граничными условиями:

на входе (г = 1) : / = 1, Н = 0, б = 0; Р = 0;

на стенках канала: (г = к ): / = 0, Н = 0, б = 1;

(ф = а): / = 0, Н = 0, б = 1;

(ф = р): / = 0, Н = 0, б = 1;

д/ дб дР

в плоскости симметрии: (г = 0) : Н = 0, — = 0,------------= 0, — = 0;

дг дг дг

на выходе (г = Я )

через отверстие: / = /, Н = 0, б = 1; на жесткой стенке: / = 0, Н = 0, б = 1.

— Я — «ср

Здесь Я = —, / = —

г0 « 0

(16)

безразмерные константы; к =

к(гог )

Р0

безразмерная функция; I =

число закрутки; Ие =

«0 Р0

число

«0 v

Рейнольдса.

Область интегрирования в новых переменных будет выглядеть следующим образом:

й = {(г,ф,г):1 <г <Я, а<ф<в, |г|<к },

причем на границе г = Я область отверстий описывается неравенствами © Проблемы энергетики, 2007, № 3-4

2 к (Я) (п ф-у (к - 1)) у (к - 1) у (к - 1/2 )

------------------------, при -----------< ф <--------------

т<<! -Л ) (п 7) п , к = 1,2...,п,(17)

2 к (Я )(к у - пф) у (к - 1/2 ) у к

------------------, при ---------------< ф <-----

у п п

а область жесткой стенки вокруг отверстий - неравенствами:

2 к (Я)(пф - у (к - 1)) у (к - 1) у (к - 1/2 )

------------------------, при -----------<ф<----------------

у п п

|г| > •{ -----.. ч .. . ч , к = 1,2..., п. (18)

2к (Я)(ку - пф) у (к - 1/2) ук

------------------, при ---------------<ф<-------

у п п

Таким образом, система уравнений (12)-(15) с соответствующими граничными условиями (16) представляет собой полную математическую формулировку задачи течения вязкой несжимаемой жидкости во вращающемся с постоянной угловой скоростью криволинейном конвергентном канале.

Для решения задачи применим метод конечных элементов на основе метода Галеркина.

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

Согласно основной идее метода решение системы (12)-(14) будем искать в виде линейной комбинации базисных функций:

/ (ф, Г, т )=^ (ф, Г, т ), а (ф,г, т )=^ (ф,г, т),

I I

Н (ф,Г, т )=£ Н1м1 (ф,Г, т)

I

или в матричной записи:

/ (ф, г, т )=[ (ф ,г, т )]{/},

а (ф гт )=[ (ф ,г,т )]}, (19)

н (ф, г , т )=[ (ф, г , т )]{н },

где [N1=[ N е (ф, г , т), N 2 (ф, г , т),..., мег (ф, г , т) ],

f1 GX H1

{f }=• f2 •, {G}=' G 2 •, {h}= H 2

Jr. Gr. Hr.

Нижние индексы означают локальную нумерацию узлов, верхний индекс e - номер элемента, r - число узлов элемента. Каждая из функций

Ne (, r, z ) принимает значение, равное 1 в вершине с номером i, и значение,

равное 0 в других вершинах. Кроме того, функция N‘e (ф,r,z ) может быть

отлична от нуля только в тех элементах, которые содержат узел i. В остальной части области й она считается равной нулю.

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

ф = ф m, Ф = Ф m + АФ , r = Гт, r = Гт + Ar ,

z = z 1 (r ), z = z2 (r ) .

Применим линейное отображение: k = 2(Ф-Фm ^Дф- 1,

П = 2(r - rm )/Д - b

z = 2(z - zi (r))/(z2 (r)- zi (r))-1.

Тогда отображенный элемент будет единичным кубом:

-1 < k < 1, -1 < П < 1, -1 < Z < 1.

Используя обозначения: kp, Лр, Zp - координаты произвольного узла, k о = kkp , П о = ППр, Z о = ZZp, запишем функции формы в локальной системе координат [5]:

а) линейный элемент

N р =1 (1+k 0 )(1 + П 0 )1 + Z 0 );

Р 8

б) квадратичный элемент:

для угловых узлов

1

Nр = — (1 + Z0)(1 + n0)(1 + k0)(Z0 + n0 + k0 - 2);

8

для узлов на ребрах 5р = о, п = ±1, С = ±1,

1 2 N в =-(1 + С о)(1 + П о)(1 - 5 2). 4

Соответствующие выражения для узлов на других ребрах получаются перестановкой координат.

Введение локальной системы координат связано со спецификой численной реализации метода конечных элементов и позволяет использовать квадратурные формулы Гаусса при вычислении коэффициентов матрицы жесткости. Рассмотрим интегральную матрицу вида

/[Р №,

V

где [Р ] зависит от Nр (г, ф, т) или их производных по глобальным координатам

Э*р Э*р Э*р

-----,------ или------.

дг Эф дт

Согласно правилу дифференцирования сложной функции имеем:

дN р дN р дф дN р дг дN р дТ

Э5 Эф Э5 дг Э5 дт 05,

дN р ЭN р дф ЭN р дг ЭN р дт

— = ----------+ -------+ —=----,

Эл дф дп дг дп д т дп

ЭN р ЭN р дф Э* р дг ЭN р дТ

— = ----------+ -------+--- ---.

дф 0£ дг Э^ дт Э£

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

1 р д 1 I р & д 1 Эф дг дт'

дф 05 Э5 Э5

ЭN р = 3 -1 0* р , 3 = Эф дг эт

дг Эп 0П дп Эп

э* р 0* р Эф Эг эт

. ЭТ . . ЭС . .ЭС ЭС ЭС _

при условии, что матрица Якоби преобразования J не вырождена.

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

йт йф йг = det (J )й5 йц .

Исходная интегральная матрица преобразуется к виду 1 1 1

I / |[(5, п, С)]det (I № йп йС ,

-1 -1 -1

которая интегрируется с помощью квадратурной формулы Гаусса [7]:

1 1 1

* I 1

-1 -1 -1

п п п

I I 11 (5, П, С )й5 йп йС = £ £ £ Wijkf (, п р Ск ).

I=0 j=0 k=0

(20)

Функции f (51, п j, С k ) вычисляются таким образом, чтобы достигнуть

наилучшей для данного количества точек точности интегрирования. Координаты гауссовских точек приведены, например, в работе [6].

Применяя метод Галеркина и производя ортогонализацию невязки к базисным функциям, получим систему уравнений

|[^ ]Т

О

X 1

f----

т т Ие

д f д f д f

— + 1О— + Н— -дт дф Ш

( 1 ^д2 f 1 д2 f 1 д2 f ^

+

X -V Кеу

дт 2 т 2Ие дф2 Ие дг

2

2

У

1 ^ /"

-|[^]Г [ X +-----------------21т 2зт ф — йО -1 /[« ]

О

- X /[* ]

О

(д 2 НЧ

О

дтдг

дфдт т Ие 5ф йО = 12 /[^]Т (т созф)йО;

д2 О 2 дО , ч_

X----------------------+1 (О + 2соз ф)тО

(21)

О

О

т дг дтдг

йО + XI /[^ ]т

О

2

д 2 О

дфдг

йО + — /[^]т

Ие О

(д2 Н 1 д2 НЧ

- + -

дт2 т2 дф 2

йО +

+

О

т Ие

- f

дН

дН дН

- 1О--------------Н-----------+

дт 5ф дг

( 11 д 2 Н'

X + — Я2,1 йО

V Ие

(22)

/[* ]

О

2

X +

V Ие У

1 д f

1т 2 дф

+ 2 f созф +

&

Ие

дО

дО дО

+ 1От — + тН-дт дф дг

+

1

3

-/[" ]т

п

- 7 Л" ]

1 п

- + -

г г 2 Ие у

З2 С Г З 2 С Г З 2 С

- + -

+

г дї дф

дф 2 Ие дг 2 Ие дї 2

йП = 7/["]т (Гзіпф)йП;

- 2Є(/ - 7зіпф)

й П -

п

Так как базисные функции в (19) нулевого порядка непрерывности (непрерывны сами функции, но не их первые производные), то в уравнениях (21)-(23) должны присутствовать производные порядка не выше первого.

Применяя формулу Остроградского-Гаусса и учитывая граничные условия (16), уравнения (21)-(23) перепишем в виде:

/[" ]

п

к 1

/—

г г Ие у

— + 7Є---------------

дГ дф

2

к +-----27Г зіпф

Ие .

/ '

Г 2 у

йП+|

п

к------

Ие

' д["]т зт

У зг дг

й П+

+ {[" ]ТН — й П + |

П дї П

1 д["]т д/ 1 д["]т д/

г 2Ие дф дф Ие дї дї

Рд["]т дН й П + к---------------------й П +

П

дї дг

П

З["т дС дф дГ

П

( 2 дС Г Ие дф

+ 7 (С + 2соз фі)гЄ

йП = 72 /["]т(Гсоеф)йП ; (24)

П

гг ,т 1 д/ ,д["]т д/ гд["]т дС

к Г["]т---------й П - ------йП - к7 ----------йП +

П Г дї П дг дг П дг дф

+

П

1 ^ЗН

_— / г Ие У дг

дН то ( 1 1 гд[" ]т дН

^ -Н й Д- к + — й П

Зф дї у V Ие) Д ® дї

- И П

^д[" т дН 1 д[" т дНЛ

+

дг дг г 2 дф дф

й П = 0;

• Зі" т д/

П

/------------й П -/["|

П Зф дг

/[" ]т (( /г--

' 2

к + —

V Иеу

1 д/

ІГ дф

+ 2 / созф

к г 1 д["1 дН й П — Г---------------й П +

7 ПГ дф дї

(25)

+

П

Ие

дС

ЗС

дС

+ 7СГ------+ ГН------+ 2С(/ - 7зіпф)

дг дф дї

-Г д["1 ЗС

й П + Г-----------й П +

П

Ие дї дї

+

+

V Г Г 2 ИеУ

д["]т ЗЄ Г д["]т дв

+

дф дф Ие дг дг

й П = 7/[" ]т (г зіп ф)й П; (26)

П

к

1

к

1

Так как исходная система уравнений нелинейна, то алгоритм ее решения будет носить итерационный характер. При этом он строится так, что на каждой итерации сомножители в конвективных членах берутся с предыдущей итерации. Подставляя (19) в уравнения (21) - (23), получим:

П

к 1 'і

/----

V гг Ие,

З["] д["] нЗ["]'( й

+ 7С------+ Н------\й П{/}+

( 1

ЗГ

дф

дї

к—

д["]т З["] ЗГ ЗГ

й П {/}+

П

1 д["]т З["] 1 д["]т д["]

Ие дф дф Ие дї дї

( 1 М"1

йП {/}- Г["]тI к+------------2її2 зіпф — йП {/}+

■’ V Ие у Г2

П

рд["]т З["] , ^ ,

+ к---------—йП {С}-7/["]

П дф ЗГ П

гд["]т д["] 2. т

+ к----------VійП {Н}= 72 /["]т (Гсозф)йП.

2 д["] Г Ие дф

+ 7 ( + 2созф)г [" ]

й П-(С}+

П

дї дг

П

(27)

к Г["]т-йП {/}-к Гд["] йП {/}-к7 ГЗ["] йП {ї}+

* г ж * ж х » — ~

П

/[" ]т

г дї (( 1

П

дї ЗГ

П

дї дф

П

г Ие

- /

д[" ] д[" ] д["]

ЗГ

-- 7С

дф

- Н-

дї

. , ( 1 |Гд["]т З["]

йП {Н}- к+—1/^------йП {Н}-

V КеУП дї дї

- И П

д["]т д["] 1 д["]т д["]

-----------------1----------------------

дГ ЗГ Г2 дф дф

((

йП {Н} = 0;

(28)

П дф дг

(

П

2 " 1 д[[" ]

Ие у ІГ2 дф

+ 2[" ]соз ф

й П {/} +

+

3

/Г-----------

Ие у

П

к — + ■

Г Г 2 Ие

дг

т

+ ЮГ^-^- + ГН+ 2["](/ - 7зіп ф)

дф

дї

1 ^З["] З["] Г З["] З["] Г д["] д["]

+-------------+

дф дф Ие дг Зг Ие Зг

дї

к .1 д["]т д["] 1Г, ч

-Г----------—---—й П {Н }= 7 Г[" ]т (Г зіп ф)й П.

7 ПГ дф дї П

й П { } + й П {С}-(29)

Здесь индексы номера итерации и номера элемента опущены. Интегрирование проводится по области элементов, так как в остальной области "р = 0.

Эту систему можно представить в виде

' [ И, А + [ 1^ , {С } + 13 . { Н II

< [ 21 | {А}+ [ 22 , {С } + [ 23 {Н } = [ 2

[ 31 . ■{А}+ [ 32 . {С } + [ 33 . {Н } = В 3.

где введены обозначения

а

к 1

/—

г г Яе у

д[л] д[л] д[л]

/

дг

+ + Н-

дф

дг

+

а

1 д[л]Г д[л] 1 д[л]Г д[л]

Яе дф дф Яе дг дг

й а

-/[«]Г

а

й а+

1

1

к-V Яе у

д[л]Г д[л]

дг дг

й а+

2

к +------2Гг з1пф

Яе

[л ]

й а;

д[л]Г д[л]

_ й а а дФ дг а

— IДл ]Г

2 д[л ]

г Ие дф

+ I (С + 2соз ф)г [Л ]

й а;

а

а

г дг

а

дг дг

[ 22 ]=-И [

а

д[Л]Г д[Л]

дг дф

й а;

а

V г Яе

- А

д[Л ] д[Л] д[Л ]

- ГС^^ - Н-дг дф дг

й а -

1 ^рд[Л]Г д[Л]

1+-V-

V ке у а

дг дг

й а -

д[л]Г д[Л] 1 д[Л]Г д[Л]

----1---------~ + ~2---------------

дг дг г2 дф дф

д[Л]Г д[Л] '

а

дф дг

й а-|[л ]Г

а

с

й а;

2

к + — V Не

1 д[Л ]

Гг дф

+ 2[Л ]соз ф

йа;

[32 ]= ПЛ] + /

' 3 -\а[л] д[N] д[N] )

А--------------+ Юг----------+ гН-------+ 2[Л](( - Г«1п ф)

Яе У дг дф дг ,

й а +

а

к 1 л ---I-------

г г 2 Яе У

д[Л]Г д[Л] г д[Л]Г д[Л] г д[л]Г д[л] +-----------------------------+

дф дф Яе дг дг Яе дг дг

й а;

2

г

г

1

[к 33 ]=- - } I п

- Г1 д[*]Т д№]

а й;

йг дФ &

В ]= 12 |[^]Т (г ео«ф)а Й; [в2 ]= I]Т (г «іпф)аЙ; [в3 ]= 0.

й й

Обозначим через {Ф} вектор неизвестных, составленный из векторов {А },{С},{Н }:

{Ф}= {/1,/2 ,...,/г &1,^ ,...,Ог ,Ні,Н2 ,...,Нг },

через [ке ] - матрицу, составленную из подматриц [кЄ ], через {ве } - вектор

правых частей.

Запишем (30) в матричном виде

[Ке ] {Ф}Т ={ве}.

Матрица [ке ] и вектор {Ве } называются локальной матрицей жесткости и локальной правой частью элемента с номером е. Далее необходимо провести суммирование локальных матриц и локальных правых частей по всем конечным элементам области. При этом результирующая система алгебраических уравнений (называемая глобальной)

[к ]{Ф}={В} (31)

будет иметь размерность Ь х Ь по числу узлов. Суммирование выполняется по следующему правилу: локальные матрицы и локальные правые части

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

Искомые функции должны удовлетворять краевым условиям (16). Естественный способ обеспечить выполнение данного требования заключается в следующем изменении глобальной системы уравнений.

Пусть для некоторого у -го узла необходимо обеспечить равенство

ф у = А. (32)

Обнулим все коэффициенты у -й строки матрицы [ к ] за исключением диагонального к у, который положим равным единице. Одновременно присвоим

коэффициенту Ьу вектора правой части значение А. Тогда измененное у -е

уравнение системы (31) будет равносильно равенству (32). Все остальные уравнения системы (31) преобразуются путем вычитания произведения ку Ф у

из Ьу и подстановкой ку = 0,г = 1,2,...,Ь,г Ф у .

Решение полученной системы алгебраических уравнений (31) позволит найти значения компонент скоростей в узлах элемента, а по формулам (19) в любой точке элемента.

Для вычисления давления применим следующую процедуру.

С помощью уравнения неразрывности со штрафным параметром (15) вычислим давление в центре каждого элемента Рс и предположим, что это давление является постоянным внутри каждого элемента и терпит разрыв первого рода на границах элементов.

Рассчитаем значения давления в узловых точках. Для этого запишем аппроксимацию давления на конечном элементе

Р=Е = [Я ]{Р}. (33)

г

Нижние индексы обозначают локальные номера узлов, {р} - означает вектор узловых значений давления, который нужно определить.

Так как давление внутри элемента считается постоянной величиной, то левую часть выражения (33) можно заменить значением давления в центре данного элемента, т.е. Рс = [V ]{Р }.

Применим метод Галеркина

/[V]Т (Рс-[V]{Р})П= 0.

о

Перепишем это уравнение в следующем виде:

/[V ]Т Рсй й = /[V ]Т [V ]{р }а Й. й й

Обозначив [ке ]= /[V]Т [V]йЙе , [ве ]= /[V]ТРсйЙе ,

Й е Й е

запишем полученную систему уравнений для элемента в виде

[ке ]{Р }=[ве ].

Производя ансамблирование локальных матриц жесткости [ке ] в

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

формирования матриц учитываются граничные условия для давления аналогично тому, как это делалось для компонент скоростей. Решая полученную систему уравнений, найдем значения давления в узлах элемента, а по формулам (33) - в любой точке элемента.

Предложенный алгоритм позволяет определить поля скоростей и давления в проточной части радиально вращающегося криволинейного конвергентного канала в зависимости от чисел закрутки и Рейнольдса.

Summary

In work the algorithm of the decision of a three-dimensional problem of current of a viscous incompressible liquid in curvilinear convergent channel is described by a method of final elements

Литература

1. Багоутдинова А.Г., Золотоносов Я.Д. Гидродинамика течения вязкой жидкости в радиально вращающемся конвергентном канале, сочлененном с кольцевой насадкой / ВИНИТИ. - М., 2007. - 65 с. - Деп. в ВИНИТИ; №72-В2007 от 29.01.2007.

2. Соколов В.И. Центрифугирование. - М.: Химия, 1976. - 408 с.

3. Лойцанский Л. Г. Механика жидкости и газа. - М.: Наука, 1978. - 736 с.

4. Пейре Р., Тейлор Т. Д. Вычислительные методы в задачах механики жидкости: пер. с англ., Ленинград.: Гидрометеоиздат, 1986. - 352 с.

5. Сегерлинд Л. Применение метода конечных элементов - М.: Мир, 1979. - 392 с.

6. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975. -

344 с.

7. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. - М.: Наука, 1987. - 600 с.

Поступила 23.01.2007