close

Вход

Забыли?

вход по аккаунту

?

О свойствах аппроксимации одной дискретной модели несжимаемой жидкости.

код для вставкиСкачать
Вычислительные технологии
Том 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
дискретного поля скорости частиц в некоторой окрестности каждого узла прямоугольной
сетки. А именно, поле скорости аппроксимируется линейной функцией
φα =
1X
mk (aα xk + bα yk + cα zk + dα − uk )2 → min .
2 k
(1)
Работа выполнена в рамках Программы интеграционных фундаментальных исследований СО РАН
(проект №1).
c Е. В. Овчинникова, А. М. Франк, 2001.
°
∗
51
Е. В. Овчинникова, А. М. Франк
52
Здесь и далее mk — постоянные массы частиц, определяемые при начальной дискретизации, xk = (xk , yk , zk ), uk — координаты и скорости частиц. Суммирование производится по
всем частицам, принадлежащим окрестности узла сетки с мультииндексом α. Под окрестостью узла Uα или частицы Uk будем в дальнейшем понимать куб с ребром 2h и центром
в узле или частице соответственно. Тогда условие соленоидальности имеет вид
a1α + b2α + c3α = 0
(2)
(верхний индекс — это номер координаты).
Подставляя решение (1) в (2), несложно получить уравнения несжимаемости в каждом
узле сетки
X
mk Dαk uk = 0,
(3)
k
представляющие собой неголономные связи, где Dαk являются явно вычисляемыми функциями координат частиц.
Для простоты будем предполагать отсутствие внешних сил. Тогда движение частиц
описывается уравнениями Лагранжа
duk X
=
Dαk λα .
dt
α
(4)
Уравнения (3), (4) совместно с уравнениями переноса частиц
(5)
ẋk = uk
и представляют дискретную модель жидкости.
Уравнения динамики однородной идеальной несжимаемой жидкости единичной плотности, аппроксимацию которых мы будем исследовать, имеют вид
div u = 0,
(6)
u̇ = − grad p.
(7)
1. Аппроксимация оператора div
Введем характерный объем частицы µ = max mk (плотность равна единице). Для упроk
щения вычислений заменим форму записи (1) на эквивалентную ей
φα =
´2
³
1 X
mk aα (xk − xα ) + bα (yk − y α ) + cα (zk − z α ) + dα − (uk − uα ) → min,
2
k
(8)
xk ∈Uα
где u = (u, v, w). Под xα подразумевается
P
Uα
m k xk /
P
mk . Аналогичные выражения спра-
Uα
ведливы для y α , z α и uα . Предполагается, что µ << h3 .
СВОЙСТВА АППРОКСИМАЦИИ ДИСКРЕТНОЙ МОДЕЛИ
53
Задача (8) сводится к решению экивалентной ей системы линейных уравнений
´
X ³
mk aα (xk − xα ) + bα (yk − y α ) + cα (zk − z α ) + dα − (uk − uα ) (xk − xα ) = 0,
Uα
³
´
aα (xk − xα ) + bα (yk − y α ) + cα (zk − z α ) + dα − (uk − uα ) (yk − y α ) = 0,
X
mk
X
³
´
mk aα (xk − xα ) + bα (yk − y α ) + cα (zk − z α ) + dα − (uk − uα ) (zk − z α ) = 0,
Uα
Uα
dα
X
mk =
Uα
X
Uα
(9)
mk (uk − uα ) = 0.
Поскольку dα = 0, определитель системы из первых трех уравнений обозначим через ∆α ,
а также введем обозначения:
X
X
2
yy
axx
=
m
(x
−
x
)
,
a
=
mk (yk − y α )2 ,
k
k
α
α
α
Uα
X
azz
α =
Uα
axz
α =
X
Uα
Uα
mk (zk − z α )2 ,
axy
α =
X
Uα
mk (xk − xα )(zk − z α ) ,
mk (xk − xα )(yk − y α ) ,
ayz
α =
X
Uα
mk (yk − y α )(zk − z α ) .
yy
zz
xy
xz
yz
Обозначим также через Axx
α , Aα , Aα , Aα , Aα , Aα
дополнения.
Используя формулу Крамера, легко убедиться, что
¯ P
¯
m (x − xα )(uk − uα )
¯ Uα k k
¯
¯
¯ P
1
¯
mk (yk − y α )(uk − uα )
a1α =
¯
∆α ¯ U α
¯
¯ P
¯
mk (zk − z α )(uk − uα )
¯
Uα
Если обозначить вектор

1
∆α
соответствующие алгебраические
¯
axy
axz
α
α ¯¯
¯
¯
¯
yy
yz ¯
aα aα ¯ .
¯
¯
¯
¯
ayz
azz
α
α ¯
xy
xz
Axx
α (xk − xα ) + Aα (yk − y α ) + Aα (zk − z α )

 xy
yz
 Aα (xk − xα ) + Ayy
α (yk − y α ) + Aα (zk − z α )


yz
zz
Axz
α (xk − xα ) + Aα (yk − y α ) + Aα (zk − z α )
через Dαk , то нетрудно видеть, что
a1α =
X
k
xk ∈Uα
mk D1αk (uk − uα ).






Е. В. Овчинникова, А. М. Франк
54
А поскольку
X
mk D1αk
k
xk ∈Uα
¯ P
¯
axz
m (x − xα ) axy
α
α
¯ Uα k k
¯
¯
uα ¯¯ P mk (yk − y ) ayy ayz
α
α
α
uα =
¯
∆α ¯ U α
¯
¯ P
¯
azz
mk (zk − z α ) ayz
α
α
¯
Uα
то
a1α =
X
¯
¯
¯
¯
¯ 0
¯
¯
¯
¯
¯
uα ¯¯
¯
0
¯=
¯ ∆α ¯¯
¯
¯
¯
¯ 0
¯
¯
¯
¯
axy
axz
α
α ¯
¯
¯
yy
yz ¯
aα aα ¯ = 0,
¯
¯
¯
ayz
azz
α
α
mk D1αk uk .
k
xk ∈Uα
Очевидно, что аналогичные равенства выполняются и для коэффициентов b2α и c3α . Определим дискретный оператор дивергенции:
X
divh u|x=xα =
mk Dαk uk .
k
xk ∈Uα
Покажем теперь, что для u ∈ C 2 (R3 ) дискретный оператор divh аппроксимирует оператор
√
div с точностью O( 3 µ) + O(h2 ). Нетрудно показать, что определитель ∆α есть O(h15 ) .
А также
X
√
mk (xk − xi )3 = O(h7 ) + O( 3 µh5 ) ,
Ui
X
Ui
√
mk (xk − xi )2 (yk − y i ) = O(h9 ) + O( 3 µh5 ) ,
X
Ui
√
mk (xk − xi )(yk − y i ) = O(h7 ) + O( 3 µh4 ) ,
X
Ui
√
mk (xk − xi )2 = O(h5 ) + O( 3 µh5 ) .
Далее, используя разложение функции u(x) в ряд Тейлора в окрестности точки xα ,
получим
¯ P
¯
¯
¯
mk (xk − xα ) (u(xα ) + (xk − xα ) · ∇ u(xα ) + ... − uα ) axy
axz
α
α
¯ Uα
¯
¯
¯
¯
¯
¯
¯
P
yy
yz
1 ¯
¯
1
m
(y
−
y
)
(u(x
)
+
(x
−
x
)
·
∇
u(x
)
+
...
−
u
)
a
a
k
k
α
k
α
α
α
α
α
α
aα =
¯ U
¯=
α
¯
¯
∆α
¯
¯
¯ P
¯
yz
zz ¯
¯
m
(z
−
z
)
(u(x
)
+
(x
−
x
)
·
∇
u(x
)
+
...
−
u
)
a
a
k
k
α
α
k
α
α
α
α
α ¯
¯
Uα
√
= u0x (xi ) + O( 3 µ) + O(h2 ).
Аналогично получается
√
b2α = vy0 (xi ) + O( 3 µ) + O(h2 ),
√
c3α = wz0 (xi ) + O( 3 µ) + O(h2 ).
СВОЙСТВА АППРОКСИМАЦИИ ДИСКРЕТНОЙ МОДЕЛИ
55
Таким образом, доказано, что
√
divh u|x=xα = div u|x=xα + O( 3 µ) + O(h2 ).
√
В дальнейшем будет показано, что разность между xα и xα есть O( 3 µ). Следователь√
но, для u ∈ C 2 разность между div u|x=xα и div u|x=xα тоже равна O( 3 µ), т. е. дискретный
оператор дивергенции можно считать определенным в узле xα с тем же порядком аппроксимации.
2. Aппроксимация оператора grad
Очевидно, что вопрос об аппроксимации уравнения (7) сводится к исследованию аппроксимации градиента давления в правой части (4). Положим
gradh p|x=xk = −
X
Dαk pα h3 .
(10)
α
xα ∈Uk
Здесь суммирование ведется по восьми угловым узлам ячейки, содержащей k-ю частицу. Поскольку в правую часть (10) помимо самой функции p входят еще неявно координаты и массы частиц, то аппроксимацию необходимо проверять на гладких решениях системы (6), (7). Ясно, что при эволюции достаточно малых жидких объемов в силу гладких
√
точных решений их диаметр будет оставаться того же порядка 3 µ, что и для начальной
дискретизации при t = 0.
2.1. Слабая аппроксимация
Покажем, что оператор (10) аппроксимирует градиент давления в слабом смысле.
Пусть объем жидкости V частично ограничен твердой стенкой, а также имеет свободную границу. Причем на стенке выполнено условие непротекания un = 0, а на свободной
границе положим давление равным нулю. Тогда для p, u ∈ C 1 операторы div и grad являются сопряженными.
Z
Z
Z
(grad p, u) = grad p · u dx = p un ds − p div u dx = −(p, div u).
V
S
V
Здесь S — кусочно-гладкая поверхность, ограничивающая объем V . Если для дискретных
операторов divh и gradh скалярные произведения определить следующим образом:
(gradh p, u)h =
X
gradh p|x=xk uk mk ,
X
pα divh u|x=xα h3 ,
k
(p, divh u)h =
α
то нетрудно видеть, что они также являются сопряженными:
(gradh p, u)h = −(p, divh u)h .
Е. В. Овчинникова, А. М. Франк
56
Докажем теперь, что для p, u ∈ C 1 дискретный оператор градиента gradh p аппроксимирует в слабом смысле оператор grad p:
√
(grad p − gradh p, u) = O(h2 ) + O( 3 µ).
Для этого мы продолжим gradh p, определенный только в точках x = xk , на все пространство. Разделим объем V на непересекающиеся окрестности частиц Ωk , используя,
например, мозаику Дирихле, и положим
gradh p|x∈Ωk = gradh p|x=xk .
Аналогично поступим и с divh u:
divh u|x∈Uα = divh u|x=xα .
Заметим, что
(gradh p, u) =
X
gradh p|x=xk
k
Z
u(x) dx =
Ωk
³
X
√ ´
√
=
gradh p|x=xk uk mk + O(µ 3 µ) = (gradh p, u)h + O( 3 µ)
k
и
(p, divh u) =
X
α
=
X
α
divh u|x=xα
Z
p(x) dx =
Uα
³
´
divh u|x=xα pα h3 + O(h5 ) = (p, divh u)h + O(h2 ).
Тогда
√
(grad p − gradh p, u) = (grad p, u) − (gradh p, u) = −(p, div u) − (gradh p, u)h + O( 3 µ) =
√
√
= −(p, div u) + (p, divh u)h + O(h2 ) + O( 3 µ) = −(p, div u) + (p, divh u) + O(h2 ) + O( 3 µ) =
√
√
√
= (p, divh u − div u) + O(h2 ) + O( 3 µ)) = (p, O(h2 ) + O( 3 µ)) + O(h2 ) + O( 3 µ) =
√
= O(h2 ) + O( 3 µ).
2.2. Сильная аппроксимация
Покажем теперь, что для поточечной аппроксимации градиента давления условия несжимаемости (3) надо ставить не в узлах некоторой сетки, а в каждой частице.
Для начала рассмотрим случай, когда частиц бесконечно много. Это позволит значительно облегчить выкладки переходом от сумм к интегралам. Позже мы вернемся к
интересующему нас случаю дискретного распределения частиц.
Докажем, что
gradh p|x=xk = grad p|x=xk + O(h2 ).
Сначала найдем xi :
Rh Rh Rh
Rh Rh Rh
x dx dy dz
ξ dξ dη dζ
(xi + ξ) dξ dη dζ
Ui
−h −h −h
−h −h −h
= xi +
= xi .
xi = R
=
(2h)3
(2h)3
dx dy dz
R
Ui
СВОЙСТВА АППРОКСИМАЦИИ ДИСКРЕТНОЙ МОДЕЛИ
57
Аналогично получается
y i = yi ,
z i = zi .
Теперь найдем значения интегралов, входящих в определитель ∆i :
Z
Z
(x − xi )(y − y i ) dx dy dz = (x − xi )(y − yi ) dx dy dz = 0,
Ui
Ui
Z
(x − xi )2 dx dy dz =
Ui
Zh Zh Zh
ξ 2 dξ dη dζ = (2h)3
h2
.
3
−h −h −h
Подставляя найденные значения, вычислим определитель ∆i :
R
R
¯ R
¯ (x − xi )2 dx
x
(x
−
(x − xi ) (z − z i ) dx
i ) (y − y i ) dx
¯U
Ui
Ui
¯ Ri
R
¯ (x − x ) (y − y ) dx R (y − y )2 dx
(y − y i ) (z − z i ) dx
¯
i
i
i
∆i = ¯
Ui
Ui
Ui
¯ R
R
R
¯ (x − xi ) (z − z i ) dx
(y
−
(z − z i )2 dx
y
)
(z
−
z
)
dx
i
i
¯
Ui
Ui
Ui
¯
¯
¯
0
0
¯
µ 2 ¶3
¯
2
h
¯
h
9
;
¯ = (2h)
0
(2h)3
¯
3
3
2 ¯
h ¯
0
(2h)3
¯
3
 xx

zx

Ai (xk − xi ) + Ayx
i (yk − y i ) + Ai (zk − z i )
1  xy
1

yy
zy

Dik =
 Ai (xk − xi ) + Ai (yk − y i ) + Ai (zk − z i )  =
2
h
∆i
yz
zz
(2h)3
Axz
i (xk − xi ) + Ai (yk − y i ) + Ai (zk − z i )
3
Для бесконечного числа частиц дискретный градиент давления будет
виде
Z
¯
2
¯
¯ (2h)3 h
¯
3
¯
¯
=¯
0
¯
¯
¯
0
¯
¯
¯
¯
¯
¯
¯=
¯
¯
¯
¯

x k − xi
yk − y i  .
zk − z i
записываться в
Dik pi dxi .
gradh p|x=xk = −
Uk
Приведем вычисления только для первой координаты градиента давления, поскольку
для остальных они будут аналогичными:
grad1h p|x=xk =
3
=−
(2h)3 h2
Z
·
¸
0
0
0
(xk − xi ) pk + px (xk )(xi − xk ) + py (xk )(yi − yk ) + pz (xk )(zi − zk ) + ... dxi =
Uk
= p0x (xk ) + O(h2 ).
Таким образом, мы получили, что для p ∈ C 2 построенный нами оператор аппроксимирует градиент давления с точностью O(h2 ).
Е. В. Овчинникова, А. М. Франк
58
Будем полагать, что в любом ограниченном объеме число частиц конечно. Чтобы отличать данный случай от предыдущего, будем все используемые обозначения помечать
индексом d. Тогда xdi имеет вид
P
P
mk (xk − xi )
m k xk
Ui
Ui
d
P
= xi +
.
xi = P
mk
mk
Ui
Ui
Приблизим суммы, входящие в данное соотношение, интегралами:
Z
X
√
mk = dx + O( 3 µh2 ) = O(h3 ),
Ui
X
Ui
Ui
mk (xk − xi ) =
Z
√
√
(x − xi ) dx + O( 3 µh3 ) = O( 3 µh3 ).
Ui
√
3
Таким образом, xdi = xi + O( µ).
Теперь выразим суммы через интегралы с остаточными членами:
Z ³
Z
√ ´³
√ ´
x − xdi + O( 3 µ) y − y di + O( 3 µ) dx dy dz =
(x − xi )(y − y i ) dx dy dz =
Ui
Ui
¸
XZ ·
¡
¡
¢
¢
√
√
√
d
d
d
d
(x − xi )(y − y i ) + O( 3 µ) y − y i + O( 3 µ) x − xi dx dy dz + O( 3 µh4 ) =
=
k Ω
k
=
X
Ui
√
mk (xk − xdi )(yk − y di ) + O( 3 µh4 ).
Аналогично
Z
Ui
¯ xx d
¯ ai
¯
¯
¯
d
d
∆i = ¯¯ ayx
i
¯
¯
¯ azx d
i
(x − xi )2 dx dy dz =
X
Ui
√
mk (xk − xdi )2 + O( 3 µh4 ),
¯ ¯ xx
√
√
√
xy
d ¯
4
4
xz
4
¯
3
3
3
axz
i
¯ ¯ ai + O( µh ) ai + O( µh ) ai + O( µh )
¯ ¯
¯ ¯ yx
√
√
√
yy
yz
yy d
yz d ¯
4
4
4
¯
3
3
3
ai
ai
¯ = ¯ ai + O( µh ) ai + O( µh ) ai + O( µh )
¯ ¯
¯ ¯ zx
√
√
4
4
d
¯ a + O(√
3
3
3 µh4 )
zz d ¯
azy
azz
azy
a
i + O( µh )
i
i + O( µh )
i
i
√
= ∆i + O( 3 µh14 ),
√
d
Amn
= Amn
+ O( 3 µh9 ),
i
i

 xx d
d
d
d
xz d
Ai (xk − xdi ) + Axy
i (yk − y i ) + Ai (zk − z i )



1 
d
d
d
yy
yz
yx
d
d
d
d
Dik = d 
Ai (xk − xi ) + Ai (yk − y i ) + Ai (zk − z i ) 
=

∆i 

d
axy
i
zy d
d
d
d
d
zz d
Azx
i (xk − xi ) + Ai (yk − y i ) + Ai (zk − z i )
¯
¯
¯
¯
¯
¯=
¯
¯
¯
¯
СВОЙСТВА АППРОКСИМАЦИИ ДИСКРЕТНОЙ МОДЕЛИ

xy
xz
Axx
i (xk − xi ) + Ai (yk − y i ) + Ai (zk − z i )

 yx
1
yz
 Ai (xk − xi ) + Ayy
=
√
i (yk − y i ) + Ai (zk − z i )
∆i + O( 3 µh14 ) 

59


¶
µ√
3 µ

+O
=

h5

zy
zz
Azx
i (xk − xi ) + Ai (yk − y i ) + Ai (zk − z i )
µ√
¶
3 µ
= Dik + O
,
h5
µ√
¶¶
Z µ
Z
3 µ
d
Dik + O
Dik p(x) dx =
p(x) dx =
h5
Uk
=
Z
Uk
Uk
¶ X
¶
µ√
µ√
3 µ
3 µ
d
d
=
.
Dik p(x) dx + O
Dik p(x) dx + O
2
h2
h
U
k
Следовательно,
gradh p|x=xk = grad p|x=xk
¶
µ√
3 µ
+ O(h2 ).
+O
h2
Таким образом,
мы
показали что для p ∈ C 2 gradh p|x=xk аппроксимирует grad p|x=xk с
µ√
¶
3 µ
точностью O
+ O(h2 ).
2
h
3. Аппроксимация уравнения неразрывности
В заключение покажем, что движение материальных частиц в заданном гладком поле
скорости всегда приводит к некоторой аппроксимации уравнения неразрывности. Для наглядности рассмотрим двумерный случай.
Рассмотрим жидкий объем V (см. рисунок), состоящий из достаточно большого фиксированного количества частиц. Он ограничен ломаной γ, соединяющей соседние граничные частицы. Обозначим через ni+ 1 вектор нормали к соответствующим ребрам ломаной.
2
Пусть границы движутся в достаточно гладком поле скорости u, d = max |ni+ 1 |, объем V
2
i
мал, а h — его диаметр, тогда
X
1X
dV
ui + ui+1
=
=
(ni− 1 + ni+ 1 )ui =
ni+ 1
2
2
2
dt
2 i
2
i
Е. В. Овчинникова, А. М. Франк
60
=
Z
2
u · ndγ + O(d ) =
γ
Z
divudV + O(d2 ) = V divu + O(V h2 ) + O(d2 ).
V
Таким образом, при d << h << 1 получается аппроксимация уравнения
1 dV
− div u = 0,
V dt
которое приводится к уравнению неразрывности с помощью соотношения
V̇
ρ̇
=− ,
ρ
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 г.
Документ
Категория
Без категории
Просмотров
3
Размер файла
201 Кб
Теги
дискретное, одной, свойства, модель, жидкости, несжимаемой, аппроксимация
1/--страниц
Пожаловаться на содержимое документа