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

Вычислительные технологии

Том 6, № 4, 2001

О СВОЙСТВАХ АППРОКСИМАЦИИ ОДНОЙ ДИСКРЕТНОЙ МОДЕЛИ НЕСЖИМАЕМОЙ ЖИДКОСТИ*

Е. В. Овчинникова Красноярский государственный технический университет, Россия

А. М. Франк

Институт вычислительного моделирования СО РАН Красноярск, Россия e-mail: frank@icm.krasn.ru

The approximation of fluid dynamics equations by a discrete model which is in fact symplified version of method of particles for incompressible fluid is studied. It is shown that when the discrete conditions of incompressibility are imposed only in the nodes of a certain grid, the momentum equation is approximated in the weak sense. The strong approximation takes place when the discrete conditions of incompressibility are satisfied for each particle.

Введение

В работах [1, 2] на основе принципа Гамильтона была построена дискретная модель несжимаемой жидкости, являющаяся в некотором смысле сильно упрощенной версией метода частиц для несжимаемой жидкости [3]. С помощью этой модели численно получен известный эффект удержания шара вертикальной струей жидкости и даже неплохое совпадение периода автоколебаний с экспериментом. Тем не менее оставался открытым вопрос о связи такой модели с классическими уравнениями гидродинамики. В настоящей работе показана аппроксимация этих уравнений.

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

фа = 1 ^2 mk(aaXk + bayk + caZk + da - ufc)2 ^ min . (1)

k

* Работа выполнена в рамках Программы интеграционных фундаментальных исследований СО РАН (проект №1).

© Е. В. Овчинникова, А. М. Франк, 2001.

Здесь и далее mk — постоянные массы частиц, определяемые при начальной дискретизации, Xk = (xk, yk, Zk), Uk — координаты и скорости частиц. Суммирование производится по всем частицам, принадлежащим окрестности узла сетки с мультииндексом а. Под окре-

стостью узла Ua или частицы Uk будем в дальнейшем понимать куб с ребром 2h и центром

в узле или частице соответственно. Тогда условие соленоидальности имеет вид

аа+ьа+сСа =0 (2)

(верхний индекс — это номер координаты).

Подставляя решение (1) в (2), несложно получить уравнения несжимаемости в каждом узле сетки

J^mk Dak Uk = 0, (3)

k

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

Для простоты будем предполагать отсутствие внешних сил. Тогда движение частиц описывается уравнениями Лагранжа

^ = £ Dak Да. (4)

а

Уравнения (3), (4) совместно с уравнениями переноса частиц

Xk = Uk (5)

и представляют дискретную модель жидкости.

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

div u = 0, (6)

U = — grad p. (7)

1. Аппроксимация оператора div

Введем характерный объем частицы ^ = max mk (плотность равна единице). Для упро-

k

щения вычислений заменим форму записи (1) на эквивалентную ей

Фа = \ ^2 m-k( aa(Xk — Xa) + ba(yk — У a) + Ca(Zk — Za) + da — (Uk — Па)) ^ min, (8)

k

Xk eUa

где U = (u,v,w). Под Xa подразумевается mk Xk/ ^2 mk■ Аналогичные выражения спра-

Ua Ua

ведливы для ya ,Za и Ua . Предполагается, что ^ << h3 .

Задача (8) сводится к решению экивалентной ей системы линейных уравнений

^т*(аа(Жй - Ха) + Ьа(ук - Уа) + Са(^ - ^- (и* - Иа)) (ж* - Жа) = 0,

иа

аа(Хк - Ха) + Ьа(ук - Уа) + Са(^* - ^) + ва - (и*. - йа)) (у* - Уа) = 0,

иа

аа(х* - Ха) + Ьа(у* - Уа) + Са(^* - ^а) + ва - (й* - йа)) (^ - ^а) = 0,

иа

^а т* = ^ т*(и* - йа) = 0.

(9)

иа

иа

Поскольку ва = 0, определитель системы из первых трех уравнений обозначим через Да, а также введем обозначения:

а0Т = ^ Ш(хк - х«)2 , аГ = ^ Ш(Ук - У«)2 ,

иа иа

аа = ^ тк(^ - ^“)2 , ^ = X! тк(хк - Х«)(Ук - Уа) ,

иа и а

(х*, - Х«)(^ - га) , а^ = ^ т(у*. - У«)(^ - га) .

ахг = аа

и а

и а

Обозначим также через , Аау , А^2 , АаУ , А^ , Аа соответствующие алгебраические дополнения.

Используя формулу Крамера, легко убедиться, что

аа =

1

да

Ет*(Х* - Ха)(и* - Иа) «

иа

Е т*(У* - Уа)(и* - Иа) аУаУ

иа

Ет*(^ - ^а)(м* - Иа) аа^ а

иа

Если обозначить вектор

1

да

( АаХ(Х* - Ха) + АаУ(У* - Уа) + Аа^(^ - ^ \

АаУ(Х* - Ха) + АаУ(У* - уа) + Аа^(^ - ^

V Аа^(Х* - Ха) + Аа^(У* - Уа) + Ааг(^ - ^а) )

через , то нетрудно видеть, что

т* ^а* (и* Иа)*

А поскольку

х "л і «а

/ v mk Dak «а = A" Aa

ka

XfcSUa

Emk (xk - Ua xa) axy aa -IZ aa о y a н з a

Emk (yk- Ua ya) «ay ayz aa «a = a: о «ay ayz aa

Emk(zk - Ua z a) «az aZZ aa о «az aZZ aa

то

*a = X! mk Dak «k.

k Xfc eUa

Очевидно, что аналогичные равенства выполняются и для коэффициентов ba и ca. Определим дискретный оператор дивергенции:

divh u|X=Xa ^ ^ mk Dak uk.

k Xk eUa

Покажем теперь, что для u Є C2(R3) дискретный оператор divh аппроксимирует оператор div с точностью O(^^) + O(h2). Нетрудно показать, что определитель Aa есть O(h15). А также

У, mk(xk - Xi)3 = O(h7) + O(^h5),

Ui

J^mk(xk - Xi)2(yk - Уі) = O(h9) + O(^h5)

Ui

^mk(xk - Xi)(yk - Уі) = O(h7) + O(^h4), У mk(xk - Xi)2 = O(h5) + O(^h5).

Ui

и*

Далее, используя разложение функции и(х) в ряд Тейлора в окрестности точки Х0 получим

Ет^(х^ - Ха) (и(Ха) + (Хй - Ха) -V м(Ха) + ... - Йа) «Г

иа

aa =

1

а:

Emk(yk - Уа) («(Ха) + (Xk - Ха) -V «(Ха) + ... - «а) aay «а^

Ua

Emk(Zk - Za) («(Ха) + (Xk - Ха) -V «(Ха) + ... - «а) «а^ «а^

Ua

«X(xi) + O(^) + O(h2).

Аналогично получается

ьа=^y(xi)+o(^+o(h2) ca=wZ(xi)+o(^+o(h2).

о

Таким образом, доказано, что

аіу^ и|х=Ха = и|х=Ха + 0(л^Д) + 0(^2).

В дальнейшем будет показано, что разность между ха и ха есть 0(). Следовательно, для и Є С2 разность между и|х=Ха и и|х=ха тоже равна 0( 3Д), т. е. дискретный оператор дивергенции можно считать определенным в узле ха с тем же порядком аппроксимации.

2. Лппроксимация оператора grad

Очевидно, что вопрос об аппроксимации уравнения (7) сводится к исследованию аппроксимации градиента давления в правой части (4). Положим

gradh p|x=xk = — ^2 D«fc P« h3. (10)

a

xa

Здесь суммирование ведется по восьми угловым узлам ячейки, содержащей k-ю частицу. Поскольку в правую часть (10) помимо самой функции p входят еще неявно координаты и массы частиц, то аппроксимацию необходимо проверять на гладких решениях системы (6), (7). Ясно, что при эволюции достаточно малых жидких объемов в силу гладких точных решений их диаметр будет оставаться того же порядка , что и для начальной дискретизации при t = 0.

2.1. Слабая аппроксимация

Покажем, что оператор (10) аппроксимирует градиент давления в слабом смысле.

Пусть объем жидкости V частично ограничен твердой стенкой, а также имеет свободную границу. Причем на стенке выполнено условие непротекания = 0, а на свободной

границе положим давление равным нулю. Тогда для p, u G C1 операторы div и grad являются сопряженными.

(gradp, u) = J gradp ■ u dx = Jp un ds — Jp div u dx = —(p, div u).

V S V

Здесь S — кусочно-гладкая поверхность, ограничивающая объем V. Если для дискретных операторов divh и gradh скалярные произведения определить следующим образом:

(gradhP, u)h = ^gradh p|x=xfc ufc

mfc,

(р, &у^ и)^ = ^ Ра &у^ и|х=ха к

а

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

(gradh p, u)h = -(p, divh u)h.

Докажем теперь, что для p, u Є C1 дискретный оператор градиента gradh p аппроксимирует в слабом смысле оператор gradp:

(grad p — gradh p, u) = 0(h2) + 0(-3^).

Для этого мы продолжим gradh p, определенный только в точках x = Xk, на все пространство. Разделим объем V на непересекающиеся окрестности частиц ^, используя, например, мозаику Дирихле, и положим

gradh p|xenfc = gгadh p|x=xfc■

Аналогично поступим и с divhu:

divh u|x€Ua = divh u|x=xa ■

Заметим, что

(grad,, p, u) = £ p|x=xt / u(x) dx =

k nfc

= J^g^h p|x=xJ uk mk + 0(^^H = (gгadh P, u)h + 0(^Й

и

(p, divh u) = J^divh u|x=x^ / p(x) dx

а U

35

У &УЛ и|х=х^ра Л,3 + 0(^5^ = (р, и)Л + 0(Л,2).

а

Тогда

(grad р — gгadh р, и) = ^гаЛ р, и) — ^гаЛН р, и) = — (р, div и) — ^гаЛН р, и)Н + 0(-3^) =

= —(р, div и) + (р, и)й + 0(Л,2) + 0( 1^Д) = —(р, div и) + (р, и) + 0(Л,2) + 0(^) = = (р, и — и) + 0(Л.2) + 0( ^)) = (р, 0(Л,2) + 0( ^)) + 0(Л,2) + 0( ^) =

= 0(й2) + 0( ^).

2.2. Сильная аппроксимация

Покажем теперь, что для поточечной аппроксимации градиента давления условия несжимаемости (3) надо ставить не в узлах некоторой сетки, а в каждой частице.

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

Докажем, что

gradН р|х=Хк = grad р|х=хк + °(^2).

Сначала найдем хг.

/г /г /г /г /г /г

/ х дьх / / / (хг + £) ^п 11!^ ^п

ц —н —н —н , —н —н —н

г / ^х ^у (2Л,)3 г (2Л,)3 г"

иг

Аналогично получается

Уі = Уі

г, = г,.

Теперь найдем значения интегралов, входящих в определитель Д,:

J (х — ж,) (у — у,) ^х ^у = J (х — х,)(у — у,) ^х ^у = 0,

и и

Н Н Н

у (х — хі) ^х ^у = J J и* —Н —Н —Н

Подставляя найденные значения, вычислим определитель Ді:

Ді

/ (ж — хі)2 ^х

/ (х — Хі) (у — Уі) <іх / (х — Хі) (г — г,) <іх

и*

и*

и*

/ (х — Хі) (У — Уі) ^х / (У — Уі)2 ^х

/ (У — Уі) (г — ^х

и і и* и*

/ (ж — Хі) (г — г,) а!х / (у — Уі) (г — г,) <іх / (г — г,)2 а!х

и* и* и

3 (2 й2 У 0 0

0 й2 (2й)3 у 0

0 0 (2й)3

й!

у

(2й)9

й2

23

1

( Аг(хк— жі) + АГ(ук— Уі) + АГ(^— ^і) \

Аху (хк — Жі) + Ауу (ук — у і ) + а^У — ^і)

\ АГ (хк — х) + Ау* (ук — у) + (^ — ^) )

,й2

хк Хі

Ук — у

(2й) — \ ^ — г,

Для бесконечного числа частиц дискретный градиент давления будет записываться в виде

gradh р|х=хк = ^ У О;к р* ^х*.

ик

Приведем вычисления только для первой координаты градиента давления, поскольку для остальных они будут аналогичными:

3

(2й)3 й2

(Хк — Хі )

ик

§Га^Н Р|х=Хк =

Рк + рХ(хк )(жі — Хк) + рУ (хк )(уі — Ук) + р^ (хй)(г, — 2^) + ...

= РХ(х) + °(й2)-

^х,-

Таким образом, мы получили, что для р Є С2 построенный нами оператор аппроксимирует градиент давления с точностью 0(й2).

3

1

Будем полагать, что в любом ограниченном объеме число частиц конечно. Чтобы отличать данный случай от предыдущего, будем все используемые обозначения помечать индексом ^. Тогда х^ имеет вид

Етк Хк

_ _Ц___________

Хі V''

2^тк

V;

Е «к (Хк - Х)

Ці___________________

Ет '

V;

Приблизим суммы, входящие в данное соотношение, интегралами:

^^тк _ ^х + 0(/дЛ2)_ 0(Л3),

V;

V;

^Шк (Хк - Хі)_ (х - Хі) а!х + 0(^Л3) _ 0(^Л3).

V;

V;

Таким образом, х? _ х* + 0().

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

(х — Хі)(у — у) ^х ^у _ ^х — Х? + 0(-/д)) (у — у? + 0(/д)) ^х ^у

V;

V;

I (х — Х?)(у — у^) + 0(^Ку — у^) + 0(^УДКХ — Х^)

_ ^Шк(хк - Х^)(Ук - у") + 0(^Дл4)

V;

Аналогично

А?

у Х 1 Х )

и V;

аХХ? а ХУ? а аХ"? аі

УХ? а УУ? а У"? щ

а"Х? а "У? а а""? а

(х — х*)2 ^х ^у _ Е Шк (Хк — Х?)2 + 0( ^ДЛ4)

V;

аХХ + 0( ^дЛ4) «Xу + 0( ^дЛ4)

«Г + 0( ч/ДЛ4) «Уу + 0( ^Д Л4)

а|Х + 0( ^Д Л4) а"У + 0( ^ДЛ4)

_ А + 0(/ДО, Атп? _ Атп + 0(^ДЛ9),

«Г + 0( ^ДЛ4) а?* + 0( /ДЛ4) а"" + 0( ^Д Л4)

Б

к

А?

/ АХХ?(хк — Х?) + АХУ?(ук — у?) + АГ^к — г?) \

АУХ?(хк — Х?) + Ауу?(ук — у?) + АУ"?(гк — г?)

V А"Х?(хк — Х?) + А"У?(ук — у?) + А-^к — г?) /

1

A* + O( ^?h14)

1

/ AXx(xk - Xi) + AXy(yk - у*) + AXz(zk - zi) \

AyX(xk — xi) + Ayy (yk — yi) + AyZ (zk — zi) + O

у AfX(xk - xi) + Af(yk - Vi) + Af (zk - zi) )

I

Dik p(x) dx

J (vt + 0(^?)) P(x) dx

Dtk p(x) dx+o (f?)=^ Ddk p(x) dx+o с

3?

h2

Следовательно,

Таким образом, мы показали что для p Є C2 grad^ p|x=xfc аппроксимирует grad p|x=xk с

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

Рассмотрим жидкий объем V (см. рисунок), состоящий из достаточно большого фиксированного количества частиц. Он ограничен ломаной 7, соединяющей соседние граничные частицы. Обозначим через п^+1 вектор нормали к соответствующим ребрам ломаной. Пусть границы движутся в достаточно гладком поле скорости и, 1 = тах |п^+11, объем V мал, а к — его диаметр, тогда

З. Аппроксимация уравнения неразрывности

X,

= у и ■ п17 + 0(12) = у &уи^ + 0(12) = V&уи + О^к2) + 0(12).

7 V

Таким образом, при 1 << к << 1 получается аппроксимация уравнения

1 IV п

— —----шу и = 0,

V и

которое приводится к уравнению неразрывности с помощью соотношения

РР = _ X

Р V’

следующего из постоянства массы жидкого объема.

В сущности, этот вывод показывает, каким образом надо определять плотность дискретной среды. Поскольку классическое определение плотности с помощью предельного перехода при V ^ 0 для дискретной модели не применимо, то для корректности этого определения здесь, как и для реальных сплошных сред, необходимо, чтобы плотность не зависела от выбора контрольного объема. Точнее, должен существовать достаточно широкий диапазон размеров контрольного объема, много больших расстояния между частицами, в котором плотность не зависит от выбора конкретного объема. Очевидно, что аналогичный результат верен и в трехмерном случае. При этом кусочно-линейную поверхность объема V можно построить с использованием, скажем, сетки Делоне.

Список литературы

[1] Франк А. М. Численное моделирование удержания шара струей жидкости // Докл. АН СССР. 1999. Т. 365, №3. C. 346-349.

[2] FRANK A. M. Discrete modelling of a liquid jet suspending a ball // Russ. J. Numer. Anal. Math. Modelling. 2000. Vol. 15, No. 2. P. 145-161.

[3] Франк А. М., Огородников Е.И. Метод частиц для несжимаемой жидкости // Докл. РАН. 1992. Т. 326, №6. С. 958-962.

Поступила в редакцию 9 января 2001 г.