Перенос завихренности: вариационный подход Текст научной статьи по специальности «Физика»

УДК 532.521

ПЕРЕНОС ЗАВИХРЕННОСТИ: ВАРИАЦИОННЫЙ ПОДХОД

А.С. Толстуха

Algorithm of the calculation of the steady fluid flow with vorticity through three dimentional chanel is prsented. The velocity field is decomposed in terms of potential and vortex-induced parts. Vorticity transport calculation is based on the approximate calculation of the stream line and the time scale functions of the flow particles.

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

Рассматривается задача о стационарном течении идеальной жидкости в пространственном канале.

В пространственной области О, ограниченной границами типа Т - стенка, Г0 - граница втекания жидкости в область, Ti - граница вытекания, требуется найти вектор-функции

и = и(х) = (ui(x),u2(x),u3(x))T , й = Д(х) = (iOi (х), сг2(х), сг3(х))Т

и скалярную функцию

Р = р(х),

где х = (ад, ж2, ж3)Т £ О; и - скорость в точке х; 33 - завихренность; 33 = rot и; р - давление, удовлетворяющие уравнениям

divu = О,

(и • V) сЗ — (сЗ • V) и = О, » dui duj

i=i j=i 3

0 1998 А.С. Толстуха

E-mail: ast@univer.omsk.su Омский государственный университет

(1)

(2)

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

и • п = 0, х £ Т

u = и0, х G Г0 и • n = £Д, х £ ГД

Vp • n = (u • V) u • n, xG Г0иТиГ1.

(4)

(5)

(6) (7)

Здесь n = n(x) - внешняя к области нормаль в точке границы х,

u0 = u0(x), U\ = £Д(х) - заданные на соответствующих границах функции с

условием совместности

J и0 ■ n dS Т J Ui dS = О,

Го Г!

причем считаем и0 • п < е < 0 в каждой точке границы Г0. Для однозначности в определении давления потребуем выполнения условия

j pdx = 0. (8)

a

Первое из уравнений (1)-(3) есть уравнение сохранения массы, второе -уравнение переноса завихренности, третье определяет поле давлений; (2) и (3) являются следствием уравнения сохранения импульса. Вместе они эквивалентны уравнениям Эйлера для идеальной жидкости. Существование и единственность решения задачи протекания завихренной жидкости через замкнутую пространственную область установлены (см. [1]).

Далее полагаем, что Г0 лежит в плоскости ж3 = 0; внешняя нормаль есть п = (0,0, —1)Т и, следовательно, и03(х) > е > 0.

2. Метод решения

Для решения предлагается итерационная процедура. Начальное приближение поля скорости взять в виде:

и(°) = УД,

где ф = </>(х) - потенциал скоростей, являющийся решением задачи

Аф = 0, х £ О

v</>- п = о, хет

УД • П = -u03, X G Го уд • п = [Д, х е гд

J Дйх = 0.

a

Пусть имеется k-ое приближение поля скорости и^Д удовлетворяющее уравнению (1) и граничным условиям (4)-(6). Тогда:

1) Решить задачу определения поля давлений р(к\ удовлетворяющего ref3,

(7), (S). _

2) Найти компоненты вихря о70 на поверхности втекания. Вместе с

_ du02 du0i

^03 — щ д

UX\ UX2

из уравнения импульсов в форме Громеки-Лэмба получим

1 f , д ( {к) 110 • и0\\

^01 = ^01^03 + т,— [Ру 1 Н-------д—

и03 \ дх2 V 2 /)

1 ( д ( (к) u0- и0\\

^02 — --- U02UJ03 — —--- I рУ ' Н----- )

u03 \ dXi V 2 / J

3) Найти вектор-функцию Х(х) = (Xi(x), Х2(х), Х3(х))Т такую, что

Г (u(fc) • V)X = 0, х G О \ Х(х0) = х0, х0 G Г0

постоянную на линиях тока

I у(0) = х0, х0 е Го

В силу выбора поверхности Г0 имеем Х3(х) = 0.

Найти функцию 0(х) - время, затраченное частицей с момента вхождения в область до попадания в точку х.

Г u(fc) • V0 = 1, хеп \ 0(хо) = 0, х0 G Го

4) По функциям Х(х), 0(х) получить матрицу J преобразования координат при движении частицы вдоль линии тока.

J=(lrMSr)+u“(X)(V0)T' (12>

Действительно, рассмотрим линию тока АВ из точки х0 G Г0 в точку х и малое приращение Ах (см. рис. 1). Нужно вычислить Ау - приращение к х0 = Х(х) такое, чтобы время, затраченное частицей на путь из точки у = х0 + Ау (обозначим её литерой С) в точку х' = х + Ах (обозначим D), в первом приближении совпало со временем tAB, затраченное частицей на путь из точки х0 в точку х. Отметим на линии тока, проходящей через точки (7, D, точку Л/, лежащую на границе и равную Xq = Х(х'), а также точки С1 и D', полученные из точек А, В приращениями, ортогональными к вектору скорости в этих точках. Заметим, что время, затраченное частицей на путь из точки С1 в точку D' совпадает с tAB и, следовательно, tec = Idd’ = At. Имеем

Ау = Дх0 + Дф u(xg) + At u (xg + Дфи(хд)),

В и

где Ах0 = Xq — х0, At0 = tA'C'- В первом приближении

Ау = Дх0 + (Дф + At) u0,

где u0 = u(X(x)). Так как время нахождения частиц в области известно, в первом приближении получим:

At0 + At = 0(х') — 0(х) = V0 • Дх.

Вычисляя пределы в направлениях Xj

г™ ДуДх + ае,) _ дуг . . _ 1 0 Q „ _ гх х г щ Iim , l,z,o, ej (oij, 02j, ^3j)

ce-t 0 a OXj

получим (12).

5) Найти приближение поля завихренности по формуле (см. [2])

J СО = CTq(X).

(13)

6) Индуцированную завихренностью скорость v = v(x) найти как значение интеграла (см. [3])

1 [ <3(6 х (х-

v х =

Й7Г

CGR3

|х - С I3

7) Следующее, (к + 1)-ое приближение поля скорости есть сумма

и{к+1) = V</> + v,

где ф = </>(х) - потенциал, являющийся решением задачи

Аф = 0 х £ О

\7ф ■ n = —(v • п), х £ Т

V</> • п = —Ноз — (v • п), х G Го

V</> • n = Ui — (v • n), х G Ti

J ф chc. = 0.

a

(14)

3. Численная реализация

Для решения используется метод конечных элементов. Область О разбивается на четырехузловые элементы - тетраэдры. Вершины тетраэдров при этом образуют узлы разбиения с координатами Хд i = 1,...,Х (N - количество узлов). На этом разбиении определяются финитные базисные функции таким образом, что если ТДх) - базисная функция, связанная с узлом г, то в любом элементе она линейна и ТДхД = 1, ТДхД = 0, при j ф- г.

Метод численного решения задачи Неймана для уравнения Пуассона, с линейными базисными функциями изложен, например, в [4]. К таковым относятся: задача нахождения начального приближения (9), задача на давление пункта 1), перехода к следующему приближению поля скорости пункта 7) итерационной процедуры.

Решение задачи определения линий тока (10) ищем в виде

N

Х(х) = ]ГХгК-(х)

г = 1

методом наименьших квадратов. Пусть Х(х) - компонента вектор-функции Х(х). Требуется найти коэффициенты Ад г = 1,. . . , N обеспечивающие минимум функционалу

Ф(А1,...,Адт) = J (и • VX)2 йх a

Условие минимума квадратичного функционала

<9Ф

0,

1,...,Х

Полагая и кусочно-линейной на конечноэлементном разбиении функцией, получим

N

aijXj = 0, i = 1,..., N (15)

j=i

ад = ^2 ( -(ит • VKt)(ut • Wyr) + — ■ VKT)(ufc • WjT)

T3i,j \ keT /

V oli

где суммирование осуществляется по тетраэдрам, содержащим в качестве вершин узлы с номерами г, j; V о1т - объем тетраэдра; = и(х&), и у = ( ^ ufc)/4.

кет

Так как на границе Г0 значения Х(х) известны, система (15) редуцируется к виду

^ ^ ^ ^ ^ ^ {Ф • • • ? Х|хг' ^ Г0}.

je{i,...,jv|x^r0} щ{1,...,м|хщг0}

Решение этой системы линейных алгебраических уравнений с симметричной матрицей может быть получено численно.

Решение задачи определения временной шкалы для частиц (11) строится аналогично. Матрица системы совпадает с матрицей задачи на линии тока, несколько иная правая часть.

aijQj = bi, г £ {1,..., 1У|хг ^ Г0}

je{i,...,iV|xjgro}

bi = У^(ит • VVlT)VolT.

Т Эг

Вычисление поля завихренности - шаги 4), 5) - производится в каждом тетраэдре. Заметим, что так как функции Х(х), 0(х) получены в среднем, локальное вычисление Jy и, следовательно, Сот происходит с ошибкой, что влечёт за собой неточное выполнение условия diver = 0. Более последовательный подход заключается в требовании выполнения условия (13) в среднем, но он приводит к системе линейных алгебраических уравнений с матрицей, имеющей размеры 3N х 3N.

При расчете индуцированной завихренностью скорости практически применялась формула

VX; =

сот

хг - Хт,

4тг

{Х||хг-хт |<г}

■V oli

Хг - хт

г = 1,...,ЛГ

с некоторым г, Ху = ( ^ xfc)/4.

кет

Чтобы обеспечить выполнение условия (1) очередным приближением поля скорости, ведущее уравнение задачи (14) шага 7) нужно заменить на

Аф = —divv, х £ О.

Рис. 2

4. Пример расчета

Область - труба радиуса 1, длиной 3. Скорость на входе

и0(х) = ( ж2/2, яд/2, 1)Т , хе Г0

нормальная составляющая на выходе

Ci(x) = 1, х £ Г\.

Точное решение:

и(х) = (-ж2/2, ад/2,1)т , Ф(х) = (О, 0,1)т , xGO

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

На рис. 2 изображены результаты расчета в среднем сечении трубы: а) расчетная сетка; б) вид поля скорости в проекции на плоскость (ж!,ж2); в) линии уровней функции щ; г) линии уровней функции и2.

Таким образом, построен и реализован на основе метода конечных элементов итерационный алгоритм, решающий нелинейную задачу расчета вихревого

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

Литература

1. Антонцев С.Н., Кажихов А.В., Монахов В.Н. Краевые задачи механики неоднородных жидкостей. - Новосибирск: Наука, 1983. 316 с.

2. Овсянников Л.В. Лекции по основам газовой динамики. - М.: Наука, 1981. 368 с.

3. Бэтчелор Дж. Введение в динамику жидкости. - М.: Мир, 1973. 758 с.

4. Faizullin R.T., Tolstukha A.S. Unsteady 2 and 3-dimensional Calculations in Cascades // Unsteady Aerodynamics and Aeroelasticity of Turbomachines / Edited by Y.Tanida, M.Namba. - Amsterdam - Lausanne - New-York - Oxford - Tokyo, 1995. P.39-53.