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

УДК 532.5

МОДЕЛИРОВАНИЕ ЛАМИНАРНОГО ТЕЧЕНИЯ ВЯЗКОЙ СЖИМАЕМОЙ ЖИДКОСТИ

ПРИ МАЛЫХ СКОРОСТЯХ С. М. Аульченко, Е. И. Васильева, В. О. Каледин

MODELLING LAMINAR FLOW OF A VISCOUS COMPRESSIBLE FLUID AT SMALL SPEEDS

S. M. Aulchenko, E. I. Vasilieva, V. O. Kaledin

Рассматривается задача моделирования стационарного течения вязкой сжимаемой жидкости в ограниченной области. Частный случай несжимаемой среды получается предельным переходом в определяющих уравнениях. Разработан алгоритм расчёта поля скоростей вязкой сжимаемой жидкости. Получено численное решение двумерной задачи. Проведено сравнение с аналитическим решением задачи течения вязкой несжимаемой жидкости в плоском канале (течение Пуазейля).

The problem of modeling a stationary flow of a viscous compressed fluid in the limited area is considered. The special case of the incompressible environment turns out limiting transition in the defining equations. The algorithm of calculation of a field of speeds of a viscous compressed fluid is developed. The numerical decision of a two-dimensional problem is received. Comparison with the analytical decision of the problem of the flow of a viscous incompressible fluid in the flat channel (Poiseuille flow) is provided.

Ключевые слова: вязкая сжимаемая жидкость, поле скоростей, алгоритм, вязкая несжимаемая жидкость.

Keywords: viscous compressed fluid, field of speeds, algorithm, viscous incompressible fluid.

Введение

Описание низкоскоростного обтекания упругих тел требует решения связанных задач гидроупругости. Известны методики, основанные на конечноразностных схемах для описания течения и конечноэлементных - для определения упругих деформаций [1, 2, 3, 4]. Необходимость сопряжения разных расчётных схем вносит в расчёт трудно устранимую погрешность, поэтому целесообразна разработка моделей и алгоритмов, позволяющих использовать единую дискретную модель. В настоящей работе рассматривается приближённый подход к построению конечноэлементной модели течения жидкости. Построение корректной численной схемы для интегрирования уравнений движения сплошной среды до настоящего времени представляет собой сложную проблему. Наряду с известными сложностями в постановке граничных условий, следует отметить зависимость решения от реологии среды, что требует анализа влияния выбора определяющих уравнений на получаемые в расчёте поля скоростей и давлений. В настоящей работе рассматривается сплошная среда без связей, определяющее уравнение которой включает слагаемые со сдвиговой и объёмной вязкостью. Это позволяет использовать наиболее общую вариационную постановку задачи и упрощает выбор граничных условий. Предельный переход от рассматриваемой модельной среды к среде с бесконечно большой объёмной вязкостью даёт несжимаемую жидкость. Поэтому в рамках предлагаемой модели возможно приближённое решение уравнений движения несжимаемой среды.

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

плоском канале. Показано, что при неограниченном увеличении коэффициента объёмной вязкости полученное решение приближается к классическому течению Пуазейля.

Математическая модель

Деформация рассматриваемой среды может быть представлена в виде суммы объёмной деформации и деформации формоизменения [5]:

еЧ =е03у + , (1)

где Вц = 0, В - девиатор тензора деформации.

Уравнение для деформации формоизменения примем совпадающим с этим уравнением для вязкой ньютоновой жидкости:

т =2 в , (2)

где /и - коэффициент сдвиговой вязкости.

Остановимся на реакции модельной среды на деформацию изменения объёма. Примем два предположения: 1) деформация изменения объёма оказывает влияние только на шаровую составляющую тензора напряжений; 2) существует некоторое равновесное для данной среды давление, при котором она не изменяет своего объёма. При уменьшении давления среда начинает расширяться с постоянной скоростью, при увеличении - сжиматься, также с постоянной скоростью.

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

\ *

Р - Р Ру + Ту - Р 5у , (3)

(р - Р*) = -3£?0, Т у = 2АёЦ - ё081] )

1

вязких напряжений, ^ = _ ~ скорость объёмной

*

деформации, % - объёмная вязкость, р - физический параметр, имеющий размерность давления, а^,ё^ -тензоры напряжений и скоростей деформации.

где

- тензор

Соотношение (3) нетрудно обратить и выразить скорости деформаций через напряжения. После выполнения такой замены получим неоднородную алгебраическую систему уравнений с тремя неизвест-

ными є xx ,є yy ,є zz :

£ + з n^xx + ^ ^ ^Jyy + ^ з м l^zz axx + p ,

і-~ м yxx + [ £+~ ^Ifyy +(і-^ м ]^zz — ayy + p , (4)

З

З

^ 3 хх + з ^^УУ + ^ + р ,

ту = 3^т

После решения этой системы, мы получим определяющие соотношения для ё:: , т. е. функциональ-

и

ные зависимости є

{°ij)

IJ V IJ

i<2oxx-ayy ~azz)+-^Pxx + ayy + azz)+ 2P м

і-

бі z) + ) )

у axx + 2ayy azz) + A°xx + ayy + azz) + 2p м

є„,---------------------------3-----------------------------, (5)

і

бім z)+2 )а

t, axx ayy + 2ozz/ + 3 W^xx + ayy + °zz) + 2p м

б£м

є-

2м

lim Exx -і^-да б

lim єУУ

І^-да

2oxx -^yy - а zz

м

1 oxx - УУ + а zz

б м

1 а xx + аУУ - 2o zz

б м

(б)

aji, j +pbi -pui ■

(В)

где Ь: - проекции объёмной силы, р - плотность

жидкости, и - вектор скорости.

Подставив определяющие уравнения (3) в уравнения движения сплошной среды (8), получим уравне-

ние движения, пригодное только для частного класса сред - вязких сжимаемых жидкостей:

pu г —(і + J ^j , ji + м° г ,jj + pb i. (9)

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

puu — udivo + p(b и). (10)

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

д д t

-+ u V u .

(11)

Левая часть уравнения (10) представляет собой

(12)

производную от кинетической энергии жидкости:

л ( 2 I

- ^ d p u

p DO - ---

dt

dK

dt

Проинтегрируем (10) по объёму области:

dQ + J pu Єї)dQ — Q

-J ( 2 л

pu

д t Q 2

—J av d Г

д Q

(1З)

или

Ґ 2 I pu

dQ + J po^udQ -Q

Учитывая, что модельная среда - сжимаемая и обратимая, проанализируем её поведение при неограниченном увеличении коэффициента объёмной вязкости £ ^ да, считая, что давление р ограничено. Тогда из (5) следует:

Кш Ezz

значит, ё0 = 0, и среда приближается по свойствам к

несжимаемой вязкой жидкости.

В пределе реологическое соотношение переходит в следующее:

аи =- р8У + 3^ч ,ёкк = 0 (7)

Зде сь р имеет смысл реакции внутренней связи и

не определяется реологией.

Уравнение движения принято в виде [6]:

— I дt Q

= | qudQ. - N и + 2щ Б---------div и I +

дО ^ 3 ))

+ | р Б0dn + | р(Ьи)О.

О 0 О

Уравнение (13) выражает баланс энергии: левая часть равна производной полной механической энергии по времени, а правая - мощности заданных сил на границе области и мощности поля давлений величины р* на деформации изменения объема.

Представим искомое поле скоростей в виде линейной комбинации базисных функций с коэффициентами, равными узловым скоростям, что приводит (13) к его дискретному аналогу:

Мй + (С + £ )и = 0, (14)

Т

где М = | N pNd0. - матрица масс;

О

с = \ pNT єNdО - матрица конвективных масс;

£2

Б = Ви - матрица-столбец скоростей деформаций;

N - матрица функций формы; и - матрица-столбец

узловых скоростей; £ = |вТDBd О - матрица демп-

О

фирования; В - матрица деформаций, содержащая частные производные по пространственным координатам; д = ; NT qdГ + / В0Т р * dQ - вектор эквива-50 О 0

лентных узловых сил, В0 - матрица, связывающая скорость объемной деформации с узловыми скоро-ЛГТ ВТ ВТ

стями; N , В , В0 - транспонированные матрицы;

(7

2

О - матрица вязкости, которая имеет следующий вид:

(

О =

#+3 ,) Г 2") Г1") 0 0 0

і „) (*+4") 0 0 0

7 ') 7") (*+і ") 0 0 0

0 0 0 и 0 0

0 0 0 0 и 0

0 0 0 0 0 и

Л

Т Т

S = IВ ББа □, М = I N pNdQ.. □ О

В частности, стационарное поле скоростей получается, если в системе уравнений (14) положить и =0. Тогда имеем нелинейную систему уравнений:

(С (и) + 5 )и = 0. (15)

Таким образом, получены разрешающие уравнения для дискретного аналога задачи гидродинамики. Целью гидродинамического расчёта является нахождение полей скоростей. Плотность и вязкость, входящие в уравнения, считаются известными.

Алгоритм расчёта узловых значений скоростей состоит в следующем.

Шаг 1. Задаются начальные значения вектора скорости и = и>о, эквивалентных нагрузок 0, шаг по времени Д/ и точность решения є . Номер итерации п положим равным нулю.

Шаг 2. Вычисляются матрица демпфирования и матрица масс, которые не пересчитываются в дальнейшем:

(16)

Шаг 3. На п-й итерации вычисляем матрицу конвективных масс, которая изменяется после каждого изменения вектора скоростей:

С = ^ гг(ип )ыаО. (17)

Шаг 4. Решением системы уравнений (14) находим следующее приближение вектора узловых скоростей:

-1 -1

ип+1 =-А/М (С + Б)ип + АМ Q + ип . (18)

Шаг 5. Проверяем условие сходимости итерационного процесса для данного шага по времени:

|^п+1 ~°п\ <Е- (19)

Шаг 6. Если сходимость не достигнута, увеличиваем п на 1 и переходим к шагу 3 для этого же интервала времени, иначе положим = ип+1 - начальное

приближение для следующего интервала времени, п=0, / = / + А/ и перейдём к шагу 3 для следующего интервала времени.

Результаты расчётов

Работоспособность алгоритма проиллюстрирована на следующем примере. Пусть в области О = {(х, {): 0 < х < 200,0 < у < 100} протекает вязкая сжимаемая жидкость (рис. 1).

Рис. 1. Область определения задачи

Поставим следующие граничные условия: на правой границе зададим давление р = 0 и скорость в

средней точке итах = 0.01, а давление на левой границе принято равным перепаду давления в течении Пуазейля при указанной максимальной скорости. Параболическая зависимость скорости от расстояния до неподвижных стенок (течение Пуазейля) имеет вид [7]:

Ар

и =-----у(100 - у). (20)

2/иЬ

Задача (15) решена с использованием метода конечных элементов [8]. Исходная область была разбита на треугольники с линейной интерполяцией скоростей. Результаты численного решения задачи приведены на рисунке 2.

При увеличении коэффициента объёмной вязкости приближаемся к несжимаемой жидкости. Отклонение численного решения задачи течения вязкой сжимаемой жидкости от аналитического решения задачи течения вязкой несжимаемой жидкости в плоском канале (течение Пуазейля) при % = 500 составляет 2 %.

Рис. 2. Численное решение задачи течения вязкой жидкости в плоском канале при стремлении коэффициента объёмной вязкости к бесконечности % ^ <х : 1) % = 10 ; 2) % = 50 ;

3) % = 100 ; 4) % = 500 ; 5) течение Пуазейля

Выводы

Таким образом, предлагаемая численная схема по зволяет рассчитывать течение вязкой жидкости. Уве личение коэффициента объемной вязкости дает воз

Литература

1. Спиридонов, А. А. Решение связанной нестационарной задачи гидроупругости методом конечных элементов: автореф. дис. ... канд. тех. наук / А. А. Спиридонов. - СПб., 1992. - 15 с.

2. Аульченко, С. М. Моделирование механизма снижения сопротивления оболочек тел вращения, обтекаемых вязкой жидкостью / С. М. Аульченко, В. О. Каледин, Ю. В. Аникина // Письма в ЖТФ. - 2007. - Т. 33. -Вып. 17. - С. 83 - 88.

3. Аульченко, С. М. Вынужденные колебания оболочек тел вращения, обтекаемых вязкой жидкостью / С. М. Аульченко, В. О. Каледин, Ю. В. Шпакова // Письма в ЖТФ. - 2009. - Т. 35. - Вып. 3. - С. 33 - 39.

4. Седова, Е. А. Решение связанной задачи гидроупругости / Е. А. Седова // IX Межрегиональная научнопрактическая конференция студентов и аспирантов: сборник трудов: в 3 т. - Т. 1. НФИ КемГУ. - Новокузнецк, 2009. - С. 8 - 11.

5. Рейнер, М. Реология / М. Рейнер. - М.: Ф-М, 1965. - 224 с.

6. Мейз, Дж. Теория и задачи механики сплошных сред / Дж. Мейз. - М.: ЛКИ, 2007. - 320 с.

7. Филиппов, Н. Н. Общая физика. Гидродинамика / Н. Н. Филиппов. - М.: МАИ, 2004. - 36 с.

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

Информация об авторах:

Аульченко Сергей Михайлович - доктор физико-математических наук, доцент, главный научный сотрудник института теоретической и прикладной механики им. С. А. Христиановича СО РАН, 8-913-450-87-08, aultch@itam. nsc.ru.

Sergey M. Aulchenko - Doctor of Physics and Mathematics, Associate Professor, Chief Researcher at S. A. Khris-tanovich Institute of Theoretical and Applied Mechanics of the Siberian Branch of the RAS.

Васильева Елена Игоревна - аспирант кафедры математики и математического моделирования Новокузнецкого института (филиала) КемГУ, 8-908-948-53-54, [email protected].

Elena I. Vasilieva - post-graduate student at the Department of Mathematics and Mathematical Modeling, Novokuznetsk Institute (branch) of Kemerovo State University.

Каледин Валерий Олегович - доктор технических наук, профессор кафедры математики и математического моделирования Новокузнецкого института (филиала) КемГУ, 8-923-460-63-43, [email protected].

Valery O. Kaledin - Doctor of Technical Sciences, Professor at the Department of Mathematics and Mathematical Modeling, Novokuznetsk Institute (branch) of Kemerovo State University.

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