close

Вход

Забыли?

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

?

Численное решение нестационарных уравнений Навье-Стокса.

код для вставкиСкачать
Вычислительные технологии
Том 13, Специальный выпуск 4, 2008
Численное решение нестационарных уравнений
НавьеСтокса
К. С. Иванов
Кемеровский государственный университет, оссия
e-mail:topspin83mail.ru
In this paper we propose a method for solution of the non-stationary Navier?Stokes
equations describing a viscous incompressible flow. Our method relies on the minimal
residuals method with a multiparametric optimization based on a componentwise minimization for the residual norm of the approximate solution.
Введение
ассмотрим в области G нестационарную систему уравнений НавьеСтокса, описывающую плоское движение вязкой однородной несжимаемой жидкости. В большинстве
случаев данную систему записывают в переменных ункция тока вихрь и на каждом шаге по времени решают сначала линеаризованное уравнение переноса вихря для
? , затем уравнение Пуассона для ункции тока ? [1?. Преимущество такой постановки
задачи относительная простота реализации численного алгоритма. Однако такому
подходу присущи и существенные недостатки: во-первых, на каждом временном шаге
приходится решать два уравнения, одним из которых является уравнение Пуассона;
во-вторых, возникают проблемы, связанные с постановкой краевых условий для вихря ? . ешению этих проблем посвящено достаточно большое количество работ (см.,
например, обзор в [2?). Менее популярны методы решения системы уравнений Навье
Стокса, записанной только относительно ункции тока ? [3?. Преимуществом такого
подхода является отсутствие каких-либо существенных проблем постановки краевых
условий для ункции тока ? . Однако в этом случае на каждом дискретном временном
шаге необходимо решать системы линейных или нелинейных алгебраических уравнений
большой размерности.
1. Постановка задачи
Система уравнений НавьеСтокса, записанная в постановке ? ? ? , имеет вид
?? ?(u?) ?(v?)
+
+
= ???,
?t
?x
?y
(1)
?? = ?,
(2)
Институт вычислительных технологий Сибирского отделения оссийской академии наук, 2008.
35
К. С. Иванов
36
начальные условия
?|t=0 = ?(x, y), x, y ? G;
(3)
краевые условия
?|?G = ?1 (x, y, t), t ? [0; T ],
(4)
??
|?G = ?2 (x, y, t), t ? [0; T ].
(5)
?n
Наряду с традиционной постановкой диеренциальной задачи будем также рассматривать систему уравнений НавьеСтокса, записанную относительно только ункции тока ? :
??
??
? ??
??? ? ?y ??
?x
(6)
+
?
= ????,
?t
?x
?y
начальные условия
? |t=0 = ?(x, y), x, y ? G;
(7)
краевые условия
? |?G = ?1 (x, y, t),
t ? [0; T ],
(8)
??
| = ?2 (x, y, t), t ? [0; T ].
(9)
?n ?G
В системах уравнений (1)(5) и (6)(9) ? > 0 коэициент вязкости; ?(x, y),
?1 (x, y, t), ?2 (x, y, t) заданные ункции своих аргументов; G выпуклая односвязная
область решения; ?G гладкая граница области G. Будем считать, что задачи (1)(5)
и (6)(9) имеют единственное решение [4?.
2. Численное моделирование
Введем в области G неравномерную по t, x, y , согласованную с границей ?G сетку Gh .
Аппроксимируя задачи (1)(5) и (6)(9) на сетке Gh некоторыми разностными схемами,
получим разностные задачи:
n+1/2
?h
? ?hn
n+1/2
n+1/2
n
+ Ln+1
+ Ln+1
;
hx ?h
hy ?h = fh
? /2
(10)
n+1/2
?hn+1 ? ?h
? /2
и
n+1/2
+ Ln+1
hx ?h
n+1
+ Ln+1
= fhn+1 ;
hy ?h
(11)
?h ?hn+1 = ?hn+1 ;
(12)
?hn+1 |?G = Mh ?hn ;
(13)
?hn+1 |?G = ?n+1
h
(14)
?h ?hn+1 ? ?h ?hn
n+1
+ Ln+1
= fhn+1 ;
h ?h
?n
(15)
Kh ?hn+1 = ghn+1.
(16)
В задачах (10)(14) и (15), (16) ?h есть некоторая аппроксимация оператора Лаплаn+1
n+1
са; Ln+1
некоторые аппроксимации конвективных слагаемых в уравнеhx , Lhy и Lh
ниях движения (1) и (6); Mh какой-либо вариант условий Тома [2?, а уравнение (16)
Численное решение нестационарных уравнений НавьеСтокса
37
n+1
есть разностный аналог краевых условий (8) и (9). Заметим, что операторы Ln+1
hx , Lhy и
Ln+1
могут быть как линейными, так и нелинейными в зависимости от способа аппрокh
симации конвективных слагаемых в (1) и (6). В случае, если для аппроксимации этих
величин используются значения на верхнем временном слое, то указанные операторы
являются нелинейными, в случае же использования для аппроксимации конвективных
слагаемых значений с нижнего временного слоя данные операторы станут линейными
и мы получим линеаризованные разностные схемы.
3. Метод решения
Независимо от способа аппроксимации конвективных слагаемых разностные задачи
(10)(14) и (15), (16) для каждого дискретного момента времени можно записать как
систему алгебраических уравнений вида
(17)
A(u, u) = f,
где u, f векторы размерности m (число узлов сетки); A(u, v) = A1 (u, v) + A2 v , A2 линейный оператор, A1 билинейное отображение, обладающее следующим свойством:
A1 (?1 u(1) + ?2 u(2) , ?1 v (1) + ?2 v (2) ) = ?1 ?1 A1 (u(1) , v (1) )+
+?1 ?2 A1 (u(1) , v (2) ) + ?2 ?1 A1 (u(2) , v (1) ) + ?2 ?2 A1 (u(2) , v (2) ),
(18)
u(i) , v (i) произвольные векторы размерности m; ?i , ?i произвольные постоянные,
i = 1, 2.
Отметим, что в случае линеаризованных разностных схем на каждом шаге по времени мы имеем систему линейных алгебраических уравнений (A1 ? 0), матрица которой,
однако, зависит от временного слоя. При этом достаточно сложно установить некоторые
свойства матрицы получившейся системы линейных алгебраических уравнений (например, неособенность и знакоопределенность), которые позволили бы применять богатый
арсенал методов для ее решения [5?. В случае же нелинейной (билинейной) системы
алгебраических уравнений (A1 6= 0) данная проблема еще более обостряется. Очевидно,
что для решения системы (17) необходимо использовать такие методы, которые позволяли бы получать ее решение с использованием минимальной инормации о свойствах
оператора A.
Независимо от того, является ли система (17) линейной или нелинейной, для ее
решения построим единообразный итерационный процесс [6?
un+1/2 = un ? ?n+1 [A(un , un ) ? f ];
(19)
un+1 = un+1/2 + ?n+1 xn , n = 1, 2 . . . ,
(20)
где xn некоторый вектор размерности m; u0 произвольное начальное приближение
из области определения оператора А; ?n+1 , ?n+1 итерационные параметры.
Пусть ?n+1 квадратная матрица с m ненулевыми элементами ?kn+1
, i, j = 1...m,
ij
k произвольное целое число от 1 до m. Перепишем (20) в виде
u
n+1
=
n+1/2
yp?1
+
m
X
i=p
?kn+1
eki ,
ij
p, j = 1, 2 . . . m,
(21)
К. С. Иванов
38
n+1/2
n+1/2
где yp?1 = un+1/2 + ?kn+1
ek1 + ... + ?kn+1
ekp?1 , y0
p?1 j
ij
ненулевой ki -й компонентой.
Введем обозначение:
n+1/2
r(i)
= un+1/2 , eki вектор с одной
h
i
n+1/2
n+1/2
ki
n+1 ki
= A yi?1 + ?kn+1
e
,
y
+
?
e
? f,
i?1
ki j
ij
n+1/2
i = 1, 2 . . . m.
(22)
n+1/2
Очевидно, что r n+1 = rm
= A(un+1 , un+1) ? f и r0
= r n = A(un , un ) ? f невязки
схемы (20). Переписывая (21) относительно нормы невязки и выбирая ?kn+1
из условия
ij
n+1/2 2
[6?, можно получить
минимума ri
n+1/2 n+1/2 ri
?
ri?1 ,
(23)
i = 1, 2 . . . m, n = 1, 2 . . .
Неравенство (23) означает, что на каждом итерационном шаге норма вектора невязки не возрастает. Необходимо отметить, что в случае линейной системы уравнений
(A1 ? 0) можно показать [7?, что kr n k ? 0, при n ? ?, т. е. итерационный процесс
(19), (20) сходится при любом начальном приближении.
Приведенный алгоритм означает, что элементы матрицы ?n+1 выбираются последовательно, исходя из условия минимума соответствующей невязки. В ряде случаев
удается использовать не последовательную, а многопараметрическую оптимизацию [7?.
В этом случае
n+1/2
ui
n+1/2
ri
n+1/2
= ui?1
n+1/2
= A(ui
(24)
+ ?kn+1
zkn1 + ... + ?kn+1
zknl ,
1
l
n+1/2
, ui
) ? f,
(25)
i = 1, 2 . . . ,
и итерационные параметры ?kn+1
. . . ?kn+1
, составляющие группу, выбираются из условия
1
l
n+1/2
глобального минимума нормы вектора ri
. Для их определения на каждом итерационном шаге необходимо решить систему алгебраических уравнений (линейных или
нелинейных в зависимости от исходного оператора), размерность которой равна числу
итерационных параметров в группе. Например, если оператор A является линейным,
то система линейных алгебраических уравнений для определения итерационных пара(n+1)
(n+1)
имеет следующий вид:
метров ?k1 . . . ?kl
(Aznk1 , Aznk1 )
? (Aznk , Aznk )
2
1
?
? ...
(Aznkl , Aznk1 )
?
(Aznk1 , Aznk2 )
(Aznk2 , Aznk2 )
...
(Aznk2 , Aznkl )
?
?
?
=?
?
...
...
...
...
? ? n+1
(Aznk1 , Aznkl )
?k1
n
n
?
?
(Azk2 , Azkl ) ? ? ?kn+1
2
? Ч ? ...
...
?kn+1
(Aznkl , Aznkl )
l
n+1/2
(Aznk1 , ri
)
n+1/2
n
(Azk2 , ri
)
...
n+1/2
(Aznkl , ri
)
?
?
?=
?
?
?
?
?.
?
(26)
Численное решение нестационарных уравнений НавьеСтокса
39
Отметим, что в большинстве случаев матрица исходного оператора имеет блочнодиагональную структуру. Например, при решении разностной задачи Дирихле для
уравнения Пуассона матрица оператора имеет вид
? (n) (n)
?
(n)
a1,1 a1,2 0
0
0
a1,6 0
0
...
? (n) (n) (n)
?
(n)
0
0
a2,7 0
... ?
? a2,1 a2,2 a2,3 0
?
?
(n)
(n)
(n)
(n)
? 0
?
a
a
a
0
0
0
a
.
.
.
3,2
3,3
3,4
3,8
?
?
(n)
(n)
(n)
? 0
0
a4,3 a4,4 a4,5 0
0
0
... ?
?
?
(n)
(n)
? 0
0
0
a5,4 a5,5 0
0
0
... ?
A=?
(27)
?,
? (n)
?
(n)
(n)
? a6,1 0
0
0
0
a6,6 a6,7 0
... ?
?
?
(n)
(n)
(n)
(n)
? 0
?
a
0
0
0
a
a
a
.
.
.
7,2
7,6
7,7
7,8
?
?
(n)
(n)
(n)
? 0
?
0
a
0
0
0
a
a
.
.
.
?
?
8,3
8,7
8,8
..
..
..
..
..
..
..
..
..
.
.
.
.
.
.
.
.
.
поэтому, если итерационные параметры сгруппировать следующим образом:
n+1
n+1
n+1
n+1
{?1n+1 , ?4n+1, ?12
, ?15
, . . .}, {?2n+1 , ?5n+1, ?13
, ?16
, . . .},
n+1
n+1
n+1
n+1
n+1
{?3n+1 , ?6n+1, ?14
, ?17
, . . .}, {?7n+1 , ?10
, ?18
, ?21
, . . .},
n+1
n+1
n+1
n+1
n+1
n+1
n+1
n+1
{?8 , ?11 , ?19 , ?22 , . . .}, {?9 , ?12 , ?20 , ?23 , . . .},
(28)
то матрицы систем линейных алгебраических уравнений для их определения будут
иметь диагональный вид и решение таких систем не представляет трудностей.
Особо отметим случай, когда оптимизация проводится по всем итерационным параn+1
метрам ?1n+1 . . . ?m
одновременно. Для определения данных параметров необходимо
решить систему алгебраических уравнений, размерность которой совпадает с размерностью исходной системы. Например, в линейном случае система уравнений для опреn+1
деления ?1n+1 . . . ?m
имеет вид
?
? ? n+1 ? ?
?
(Azn1 , r n+1/2 )
(Azn1 , Azn1 ) (Azn1 , Azn2 ) . . . (Azn1 , Aznm )
?1
? (Azn2 , Azn1 ) (Azn2 , Azn2 ) . . . (Azn2 , Aznm ) ? ? ?2n+1 ? ? (Azn2 , r n+1/2 ) ?
?
?Ч?
?=?
? . (29)
? ...
? ? ...
? ? ...
?
...
... ...
n
n
n
n
n
n
n+1
n
n+1/2
(Azm , Az1 ) (Azm , Az2 ) . . . (Azm , Azm )
?m
(Azm , r
)
В общем случае решение системы представляет не меньше трудностей, чем решение
?
исходной системы уравнений. Однако, если выбрать векторы z1? , z2? . . . zm
так, чтобы
?
?
n+1
(Azi , Azj ) = 0, i 6= j , то матрица системы уравнений для определения ?1n+1 . . . ?m
будет иметь диагональный вид и определение точного решения данной системы уже не
представляет труда. Более того, можно показать [7?, что в линейном случае при таком
выборе итерационных параметров схема (19), (20) будет сходиться к точному решению
за одну итерацию.
Алгоритм групповой оптимизации по всем параметрам одновременно оказывается
весьма полезным при проведении серийных расчетов нестационарных задач в линейном случае. Например, если использовать разностную схему (10)(14) для решения
нестационарной системы уравнений НавьеСтокса, то можно заметить, что на каждом
дискретном временном шаге возникает необходимость в решении разностного уравнения Пуассона, где разностный оператор не зависит от временного шага. Поэтому для
решения данной задачи можно использовать алгоритм с полной групповой оптимизацией, затратив на первом дискретном временном шаге некоторое время на построение
К. С. Иванов
40
?
указанной выше системы векторов z1? , z2? . . . zm
. На каждом последующем временном
шаге, используя данную систему векторов, получим схему, сходящуюся за одну итерацию к точному решению. Также при проведении новых расчетов можно использовать
?
, если, конечно, при этом не
уже полученную однажды систему векторов z1? , z2? . . . zm
изменяется структура разностного оператора (например, расчеты при различных коэициентах вязкости).
Если система (17) нелинейная, то в случае плохой сходимости схемы (19), (20) для
нее аналогично линейному случаю [8? можно построить процедуру ускорения сходимости, суть которой заключается в комбинации приближений схемы (19), (20) на n-м и
(n + 2)-м итерационных шагах:
xn+2 = (1 + ?n )un+2 ? ?n un ,
(30)
где un+2, un приближения схемы (19), (20), а ?n выбирается из условия [7?
min || rn+2 || = min ||A(xn+2 , xn+2 ) ? f ||.
(31)
Заключение
Для оценки предложенного метода решения задач (1)(5) и (6)(9) проведены численные расчеты классической модельной задачи о течении вязкой однородной несжимаемой жидкости в квадратной каверне с неравномерно движущейся верхней крышкой
и задачи об обтекании вязкой однородной несжимаемой жидкостью обратного уступа.
асчеты проводились при различных значениях коэициента вязкости и при различных ункциях скорости движения крышки и входного потока. Проведенные расчеты
показали эективность предложенных алгоритмов. При достаточно больших значениях шага по времени получены устойчивые результаты. Также обнаружены нестационарные решения данных задач при стационарных краевых условиях.
Список литературы
[1?
Лойцянский Л..
Механика жидкости и газа. М.: Наука, 1987. 840 с.
[2?
оуч П.
[3?
Beam R.M.
[4?
Ладыженская А.О.
[5?
Самарский А.А., Николаев Е.С.
[6?
Захаров Ю.Н., Егорова Е.Ф., Толстых М.А., Шокин Ю.И.
[7?
Захаров Ю.Н.
[8?
Николаев Е.С.
Вычислительная гидродинамика: Пер с англ. М.: Мир, 1980. 616 с.
Newton's methods for the NavierStokes equations // Comput. Meh. 1988. Vol.
2. P. 51.II.151.II.4.
М.: Наука, 1970. 340 с.
Математические вопросы динамики вязкой несжимаемой жидкости.
вузов М.: Наука, 1978. 592 с.
Методы решения сеточных уравнений: Учебник для
Метод минимальных
невязок решения одного класса нелинейных уравнений. Красноярск, 1991 (Препр. ќ 9).
радиентные итерационные методы решения задач гидродинамики. Новосибирск: Наука, 2004. 238 с.
Нелинейное ускорение двухслойных итерационных методов вариационного типа // Журн. вычисл. математики и мат. изики. 1976. ќ 6. С. 1387.
Поступила в редакцию 20 евраля 2008 г.
Документ
Категория
Без категории
Просмотров
9
Размер файла
181 Кб
Теги
решение, уравнения, стокса, навье, нестационарные, численного
1/--страниц
Пожаловаться на содержимое документа