бусинеск. Ипм им. М. В. Келдыша ран Электронная библиотека

Название	Ипм им. М. В. Келдыша ран Электронная библиотека
Дата	19.11.2022
Размер	2.72 Mb.
Формат файла
Имя файла	бусинеск.pdf
Тип	Документы #799578

ИПМ им.М.В.Келдыша РАН • Электронная библиотека
Препринты ИПМ • Препринт № 50 за 2017 г.
ISSN 2071-2898 (Print)
ISSN 2071-2901 (Online)
Шатров О.А.,
Щерица О.В.
,
Мажорова О.С.
Опыт использования метода
Гаусса для решения разностных уравнений
Навье-Стокса в переменных
«функция тока — вихрь»
Рекомендуемая форма библиографической ссылки: Шатров О.А., Щерица О.В., Мажорова
О.С. Опыт использования метода Гаусса для решения разностных уравнений Навье-Стокса в переменных «функция тока — вихрь» // Препринты ИПМ им. М.В.Келдыша. 2017. № 50. 24 с.
doi:
10.20948/prepr-2017-50
URL:
http://library.keldysh.ru/preprint.asp?id=2017-50

Ордена Ленина
ИНСТИТУТ ПРИКЛАДНОЙ МАТЕМАТИКИ
им. М.В. Келдыша
Российской академии наук
О.А. Шатров, О.В. Щерица, О.С. Мажорова
Опыт использования метода Гаусса для решения
разностных уравнений Навье–Стокса в переменных
"функция тока – вихрь"
Москва — 2017

Шатров О.А., Щерица О.В., Мажорова О.С.
Опыт использования метода Гаусса для решения разностных уравнений
Навье–Стокса в переменных "функция тока – вихрь"
Проведён анализ эффективности использования прямых методов для ре- шения разностных уравнений, возникающих при аппроксимации уравнений
Навье-Стокса методом конечных разностей. В программном комплексе ис- пользовалась библиотека от компании Intel: Intel Math Kernel Library, ко- торая содержит параллельный алгоритм LU разложения разреженных мат- риц. В качестве теста рассматривалась задача о развитии конвекции Рэлея-
Бенара вблизи порога устойчивости. Результаты тестирования показали, что современные методы решения линейных уравнений с ленточной структурой позволяют конструировать быстрые и надёжные разностные алгоритмы для уравнений Навье-Стокса, основанные на совместном решении уравнений для функции тока и завихренности.
Ключевые слова: конвекция Рэлея–Бенара, уравнения Навье–Стокса, ма- тематическое моделирование.
Oleg Alexandrovich Shatrov, Olga Vladimirovna Shcheritsa,
Olga Semenovna Mazhorova
The experience of using the Gauss method for solving finite difference schemes of Navier-Stokes equations in the stream function–vorticity formulation
The efficiency of the use of direct methods for solving the Linear finite difference schemes for Navier—Stokes equations was studied. The software package used the Intel Math Kernel Library, which contains a parallel algorithm for LU factorization of matrices, to solve systems of linear equations with sparse matrices. The algorithms were used to solve the problem of Rayleigh-
Benard convection in a rectangular region for Rayleigh numbers near critical values. Мodern methods for solving linear equations with band structure allow to construct fast and robust algorithms for the Navier-Stokes equations based on the coupled solution of the equations for the stream function and vorticity.
Key words: Rayleigh–B´enard convection, Navier–Stokes equations, mathe- matical modeling
Работа выполнена при поддержке гранта РФФИ 15-01-03741 и Программы фундаментальных исследований Президиума РАН П-5.

3
1. Введение
Работа посвящена исследованию численных методов решения уравнений
Навье–Стокса в переменных "функция тока — вихрь". Основное внимание уделяется методам решения нестационарных задач динамики вязкой жидко- сти, когда искомые функции необходимо определить на большом промежутке времени: t >> t
ν
= H
2
/ν, где t
ν
— масштаб времени, связанный с диссипатив- ными гидродинамическими процессами, H — характерный линейный размер,
ν
— коэффициент кинематической вязкости. Примером такого процесса мо- жет служить конвекция Рэлея–Бенара — одно из самых распространенных явлений в природе и технике. Она наблюдается в атмосфере Земли и земной коре, морях и океанах, во многих технологических процессах. Характерной особенностью данного класса задач является наличие нескольких времен- ных масштабов, связанных с различными диссипативными процессами: вяз- костью, теплопроводностью, диффузией. Эти масштабы могут отличаться на несколько порядков, поэтому течение описывается жесткой системой урав- нений. Именно этим обусловлены трудности расчета длительных нестацио- нарных процессов. Используемые алгоритмы должны обеспечивать точность расчета быстрых процессов, не требуя при этом больших затрат машинно- го времени на расчет всей задачи в целом. Важно, чтобы шаг по времени определялся свойствами физического процесса, а не внутренними свойства- ми самого расчетного алгоритма.
Моделирование процессов тепломассопереноса в вязкой несжимаемой жид- кости в двумерном приближении часто осуществляется на основе уравнений
Навье–Стокса, записанных в переменных "вихрь ω, функция тока ψ" [1, 2].
Метод численного решения уравнений Навье–Стокса в переменных (ω, ψ)
можно организовать следующим образом: на каждом временном слое урав- нения для ω и ψ решаются последовательно, причем в граничном условии для ω функция тока берется с нижнего временного слоя. Такой подход обес- печивает получение сравнительно простого алгоритма, однако он неизбеж- но приводит к ограничению на величину шага по времени τ
≤
h
2
min
Re, где
h
min
= min(h
x
, h
y
), h
x
, h
y
— шаги разностной схемы по пространству, Re —
число Рейнольдса [1, 3, 4]. Это ограничение не связано со свойствами тече- ния и обусловлено только внутренними свойствами разностного алгоритма.

4
Такое ограничение оправдано при исследовании быстроменяющихся течений со сложной структурой, расчеты медленных течений, например конвекции
Рэлея–Бенара вблизи критических чисел Рэлея разумно вести с более круп- ным шагом по времени. Ограничение становится особенно обременительно,
когда решение нужно определять на отрезке времени, существенно превос- ходящем время вязкой диссипации.
Методы, в которых уравнение переноса вихря и уравнение Пуассона для функции тока решаются совместно, обладают большим запасом устойчивости и позволяют вести расчеты с шагом по времени, величина которого опреде- ляется из соображений точности решения задачи, а не требованиями устой- чивости вычислительного алгоритма [1, 5–8]. Вместе с тем использование совместного метода связано с необходимостью решать на каждом шаге по времени систему сеточных уравнений с матрицей необычной структуры. Для этого можно использовать итерационные методы [5] или прямые алгорит- мы [6]. В обоих случаях процедура решения сеточных уравнений снижает эффективность совместных алгоритмов. Хотя использование крупных шагов по времени позволяет заметно сократить время расчёта всей задачи, вычис- ление каждого шага становится дорогим как в смысле затрат оперативной памяти, так и в смысле числа арифметических операций на каждом слое.
В данной работе для решения нестационарных задач, в которых искомые функции необходимо определять на промежутке времени t >> t
ν
, применя- лись совместный и последовательный алгоритмы решения уравнений Навье–
Стокса, записанных в переменных "вихрь ω, функция тока ψ". Для решения соответствующих линейных сеточных уравнений использовалась библиотека от компании Intel: Intel Math Kernel Library, которая содержит эффектив- ный параллельный алгоритм LU разложения матриц для решения систем линейных уравнений с разреженными матрицами. Это позволило сократить затраты машинного времени на расчет одного временного слоя. Однако сов- местный метод решения уравнений Навье–Стокса все равно требует б´ольших затрат по времени и памяти для расчета одного шага по времени, чем после- довательный алгоритм. Практически безусловная устойчивость совместного метода в сочетании с современным алгоритмом решения систем линейных уравнений, учитывающим структуру матрицы системы и архитектуру вычис- лительной техники, позволяет вести расчет всей задачи значительно быстрее,

5
чем последовательный метод, в том числе на подробных сетках по простран- ству.
2. Постановка задачи
В прямоугольной области Ω = [0, L]
×
[0, H ] рассмотрим движение вязкой несжимаемой жидкости, подогреваемой снизу. Уравнения термогравитацион- ной конвекции в приближении Обербека-Буссинеска в переменных "функция тока – вихрь" можно записать в виде [9]:
∂ω
∂
t
+
∂
(ωV
x
)
∂
x
+
∂
(ωV
y
)
∂
y
= ∆ω +
Ra
Pr
∂
T
∂
x
,
(1)
−
ω
= ∆ψ,
(2)
V
x
=
∂ψ
∂
y
,
V
y
=
−
∂ψ
∂
x
,
(3)
∂
T
∂
t
+
∂
(V
x
T )
∂
x
+
∂
(V
y
T )
∂
y
=
1
Pr
∆
T ,
(4)
V
|
∂Ω
= 0,
(5)
T
|
y=0
= T
bot
,
T
|
y=H
= T
top
,
(6)
1
Pr
∂
T
∂
n
x=0,L
= 0.
(7)
V
|
t=0
= 0,
T
|
t=0
=
T
top
−
T
bot
H
y + T
bot
,
(8)
Здесь t — время, x, y — декартовы координаты,
V (t, x, y) = (V
x
,V
y
) — вектор скорости, T (t, x, y) — температура, ψ — функция тока, ω — вихрь, ∆f –
оператор Лапласа в декартовой системе координат.
Система уравнений (1) – (8) содержит два безразмерных параметра: Ra =
gβ∆T H
3
/(νκ) — число Рэлея; Pr = ν/κ — число Прандтля, где β — коэф- фициент теплового расширения, g — модуль ускорения свободного падения,
∆
T = T
bot
−
T
top
— перепад температуры, ν — коэффициент кинематической вязкости, κ — коэффициент температуропроводности.
Систему уравнений (1) – (8) необходимо дополнить граничными усло- виями для ψ и ω. На границе области заданы условия прилипания (5), с помощью соотношения (3) можно определить ψ с точностью до константы

6
и ее производную по нормали, т.е. ψ = const и
∂ψ
∂
n
= 0, обычно ψ = 0.
Физически осмысленные граничные условия для ω отсутствуют. Однако из уравнения (2) вытекает интегральное условие на завихренность:
−
ˆ
Ω
ω
dxdy =
ˆ
Ω
∆ψ
dxdy =
ˆ
∂Ω
∂ψ
∂
n
dl = 0,
(9)
которое свидетельствует о сохранении интеграла от завихренности.
Кинетическая энергия в системе изменяется по закону
∂ε
∂
t
=
−
ˆ
Ω
ω
2
dxdy,
ε
=
1 2
ˆ
Ω

V
2
x
+ V
2
y

dxdy.
(10)
Соотношения (10) получаются после умножения (1) на ψ и интегрирования полученного уравнения по области Ω. Баланс энергии (10) обеспечивается выполнением равенства
ˆ
Ω
ψ
∂ (ωV
x
)
∂
x
+
∂
(ωV
y
)
∂
y

dxdy = 0,
(11)
которое означает, что конвективные члены не влияют на баланс кинетиче- ской энергии.
Сохранение интеграла от завихренности и баланс кинетической энергии являются важнейшими свойствами течения, поэтому естественно требовать,
чтобы разностные аналоги интегральных соотношений (10), (11) выполня- лись в дискретной модели [1, 10, 11].
3. Разностная схема
Введем в области Ω прямоугольную сетку Ω
h
= Ω
h
x
×
Ω
h
y
, где Ω
h
x
=
{
x
0
=
0 < x
1
< x
2
< . . . < x
N
x
= L
}
, Ω
h
y
=
{
y
0
= 0 < y
1
< y
2
< . . . < y
N
y
= H
}
одномерные сетки по переменным x и y. Обозначим α
i+1/2
= (α
i+1
+α
i
) / 2, i =
0, 1, . . . , N
α
−
1, где α = x, y. Шаги разностной сетки:
h
α
i+1/2
= α
i+1
−
α
i
,
i = 0, 1, . . . , N
α
−
1,
h
α
i
= (h
α
i+1/2
+ h
α
i−1/2
) / 2,
i = 1, 2, . . . , N
α
−
1,
h
α
0
= h
α
1/2
/ 2,
h
α
N
α
= h
α
N
α
−1/2
/ 2,

7
где α = x, y.
Ячейки с центром в точках (x
i
, y
j
) обозначим Ω
i,j
(рис. 1-2). Функции ψ,
ω
и T заданы в узлах сетки Ω
h
. Скорости V
x
определены в точках (i, j + 1/2),
V
y
— в точках (i + 1/2, j). Сетка по времени Ω
t
=
{
t
k
= kτ, k = 0, 1, 2, . . .
}
Рис. 1. Внутренние узлы
Рис. 2. Приграничные ячейки
Разностная схема для задачи (1) – (8) получена с помощью интегро - интерполяционного метода. Проинтегрируем уравнение переноса вихря (1)
по ячейке Ω
i,j
и интервалу времени [t
k
, t
k
+ 1]:
t
k+1
ˆ
t
k
ˆ
Ω
i,j
∂ω
∂
t
dxdydt +
t
k+1
ˆ
t
k
ˆ
Ω
i,j
K
(V , ω) dxdydt =
=
t
k+1
ˆ
t
k
ˆ
Ω
i,j
∆ω
dxdydt +
t
k+1
ˆ
t
k
ˆ
Ω
i,j
Ra
Pr
∂
T
∂
x
dxdydt,
(12)
где
K
(V , ω) =
∂
∂
x
W
1, x
+
∂
∂
y
W
1, y
,
W
1, α
= ωV
α
,
α
= x, y.
Вычислим каждый интеграл в соотношении (12) по отдельности.
t
k+1
ˆ
t
k
ˆ
Ω
i,j
∂ω
∂
t
dxdydt
≈
(
b
ω
i,j
−
ω
i,j
)h
x
i
h
y
j

8
Интеграл от конвективных членов вычислим по формуле Остроградского –
Гаусса:
t
k+1
ˆ
t
k
ˆ
Ω
i,j
∂
∂
x
W
1, x
+
∂
∂
y
W
1, y

dxdydt =
t
k+1
ˆ
t
k
˛
∂Ω
i,j
−
W
1, y
dx + W
1, x
dy dt
≈
≈
τ
h
(W
1, x
i+1/2,j
−
W
1, x
i−1/2,j
)h
y
j
+ (W
1, y
i,j+1/2
−
W
1, y
i,j−1/2
)h
x
i
i
(σ
1
,σ
2
)
,
(13)
где σ
1
, σ
2
— весовые коэффициенты, σ
1
соответствует функции тока, σ
2
—
соответствует вихрю. В выражении (13) значения потоков находятся в узлах,
в которых соответствующие значения скорости и вихря не определены. Для их нахождения воспользуемся первой схемой Аракавы [1, 10]:
W
1, x
i+1/2,j
= 0.5(ω
ij
V
x,i,j
+ ω
i+1,j
V
x,i+1,j
),
W
1, y
i,j+1/2
= 0.5(ω
ij
V
y,i,j
+ ω
i,j+1
V
y,i,j+1
),
(14)
где
V
x,i,j
=
ψ
i,j+1
−
ψ
i,j−1 2h
y
j
= ψ
◦
y
,
V
y,i,j
=
−
ψ
i+1,j
−
ψ
i−1,j
2h
x
i
=
−
ψ
◦
x
Использование аппроксимации (14) обеспечивает выполнение условия (11) в дискретном случае. При аппроксимации оператора Лапласа в (12) получим:
t
k+1
ˆ
t
k
ˆ
Ω
i,j
∆ω
dxdydt
≈
τ
"
ω
i+1,j
−
ω
i,j
h
x
i+1/2
−
ω
i,j
−
ω
i−1,j
h
x
i−1/2
!
h
y
j
+
+
ω
i,j+1
−
ω
i,j
h
y
j+1/2
−
ω
i,j
−
ω
i,j−1
h
y
j−1/2
!
h
x
i
#
(σ
3
)
,
(15)
где σ
3
— весовой коэффициент. Слагаемое, описывающее подъемную силу,
аппроксимируется стандартным образом:
t
k+1
ˆ
t
k
ˆ
Ω
i,j
Ra
Pr
∂
T
∂
x
dxdydt
≈
τ
Ra
Pr
T
i+1j
−
T
i−1j
2h
x
i
h
x
i
h
y
j
Разностную схему для уравнения (1) можно записать в виде:
ω
t
+

(ψ
◦
y
ω
)
◦
x
−
(ψ
◦
x
ω
)
◦
y

(σ
1
,σ
2
)
= (ω
¯xx
+ ω
¯yy
)
(σ
3
)
+
Ra
Pr
T
◦
x

9
Для построения аппроксимации уравнения (2) проинтегрируем его по ячейке Ω
i,j
. Получим:
−
ω
(σ
3
)
i,j
h
x
i
h
y
j
=
"
ψ
i+1,j
−
ψ
i,j
h
x
i+1/2
−
ψ
i,j
−
ψ
i−1,j
h
x
i−1/2
!
h
y
j
+
+
ψ
i,j+1
−
ψ
i,j
h
y
j+1/2
−
ψ
i,j
−
ψ
i,j−1
h
y
j−1/2
!
h
x
i
#
(σ
4
)
(16)
Поделив (16) на площадь ячейки h
x
i
h
y
j
, получим:
−
ω
(σ
3
)
= (ψ
xx
+ ψ
yy
)
(σ
4
)
,
0 < i < N
x
, 0 < j < N
y
Построим аппроксимацию граничных условий для функций ψ и ω. Для этого проинтегрируем уравнение (2) по приграничной ячейке (рис. 2), для определенности рассмотрим ячейку, прилегающую к нижней границе. В ре- зультате получим граничное условие Тома [2, 12]
−
ω
(σ
5
)
i,0
h
y
0
=
"
ψ
i,1
−
ψ
i,0
h
y
1/2
#
(σ
6
)
(17)
Аналогичные условия выводятся и на остальных границах:
x = 0
:
−
ω
(σ
5
)
0,j
h
x
0
=
"
ψ
1,j
−
ψ
0,j
h
x
1/2
#
(σ
6
)
,
x = L :
ω
(σ
5
)
N
x
,j
h
x
N
x
=
"
ψ
N
x
,j
−
ψ
N
x
−1,j
h
x
N
x
−1/2
#
(σ
6
)
,
y = H :
ω
(σ
5
)
i,N
y
h
y
N
y
=
"
ψ
i,N
y
−
ψ
i,N
y
−1
h
y
N
y
−1/2
#
(σ
6
)
(18)
В угловых точках ω
0,0
= ω
0,N
y
= ω
N
x
,0
= ω
N
x
,N
y
= 0. Эти равенства полу- чаются, если проинтегрировать уравнение (2) по соответствующей угловой ячейке. Благодаря граничным условиям Тома, сеточный аналог интегрально- го условия (9) выполнен во всех узлах сетки при σ
5
= σ
6
Для построения разностной схемы для уравнения теплопереноса (4) про- интегрируем его по ячейке Ω
i,j
на интервале времени [t
k
, t
k+1
]:
t
k+1
ˆ
t
k
ˆ
Ω
i,j
∂
T
∂
t
dxdydt +
t
k+1
ˆ
t
k
ˆ
Ω
i,j
K
(V , T ) dxdydt =
t
k+1
ˆ
t
k
ˆ
Ω
i,j
∆
T dxdydt,
(19)

10
где
K
(V , T ) =
∂
∂
x
W
2,x
+
∂
∂
y
W
2,y
,
W
2,α
= V
α
T ,
α
= x, y.
Последовательно вычислим интегралы, входящие в (19).
t
k+1
ˆ
t
k
ˆ
Ω
i,j
∂
T
∂
t
dxdydt
≈
(b
T
i,j
−
T
i,j
)h
x
i
h
y
j
Конвективные члены в (19) аппроксимируются аналогично (13):
t
k+1
ˆ
t
k
ˆ
Ω
i,j
∂
∂
x
W
2,x
+
∂
∂
y
W
2,y
dxdydt
≈
≈
τ
h
(W
2,x
i+1/2,j
−
W
2,x
i−1/2,j
)h
y
j
+ (W
2,y
i,j+1/2
−
W
2,y
i,j−1/2
)h
x
i
i
(σ
7
,σ
8
)
,
где σ
7
соответствует функции тока, σ
8
— температуре. Для аппроксимации потоков W
2,x
i±1/2,j
и W
2,y
i,j±1/2
используется вторая схема Аракавы [1, 10]:
W
2,x
i+1/2,j
= V
x,i+1/2,j
T
i+1,j
+ T
i,j
2
,
W
2,y
i,j+1/2
= V
y,ij+1/2
T
i,j+1
+ T
i,j
2
,
где компоненты вектора скорости определяются следующими выражениями:
V
x,i+1/2,j
=
1 2
ψ
i+1,j+1
−
ψ
i+1,j−1 2h
y
j
+
ψ
i,j+1
−
ψ
i,j−1 2h
y
j
!
,
V
y,i,j+1/2
=
−
1 2
ψ
i+1,j+1
−
ψ
i−1,j+1 2h
x
i
+
ψ
i+1,j
−
ψ
i−1,j
2h
x
i

Аппроксимация диссипативных членов в (19) строится аналогично (15). В
итоге разностная схема для уравнения (4):
T
t
+
1
h
x
i
h
y
j
h
(W
2,x
i+1/2,j
−
W
2,x
i−1/2,j
)h
y
j
+ (W
2,y
i,j+1/2
−
W
2,y
i,j−1/2
)h
x
i
i
(σ
7
,σ
8
)
=
=
1
Pr
T
¯xx
+ T
¯yy

(σ
9
)
,
i = 1, N
x
−
1, j = 1, N
y
−
1.
Аппроксимацию граничных условий (7) можно получить, проинтегрировав
(4) по приграничной ячейке.

11
4. Аппроксимация по времени
Рассмотрим схему с весами для системы уравнений конвекции Рэлея-
Бенара:
ω
t
+
K
(ψ
(σ
1
)
, ω
(σ
2
)
) = ∆ω
(σ
3
)
+
Ra
Pr
∂
T
∂
x
,
(20)
ω
(σ
3
)
=
−
∆ψ
(σ
4
)
,
(21)
G(ψ
(σ
6
)
, ω
(σ
5
)
) = 0,
(22)
T
t
+
K
(ψ
(σ
7
)
, T
(σ
8
)
) =
1
Pr
∆
T
(σ
9
)
,
(23)
где G(ψ
(σ
6
)
, ω
(σ
5
)
) — определяет значение вихря на границе по формулам Тома
(17), (18). Система уравнений (20) – (23) является нелинейной. В данной ра- боте система (20) – (23) линеаризуется путем выбора σ
1
= σ
7
= 0. Такой под- ход позволяет избежать использования итерационного процесса на каждом временном шаге. Весовые коэффициенты σ
2
= σ
3
= σ
4
= σ
5
= σ
8
= σ
9
= 1.
Для решения системы уравнений (20) – (23) использовался метод расщеп- ления по физическим процессам [13]: сначала из уравнений Навье-Стокса
(20) – (22) определялось поле скоростей, затем вычислялось распределение температуры. В работе рассмотрены три алгоритма для уравнений Навье-
Стокса (20) – (22).
•
Метод I. Уравнения (20) – (22) решаются совместно. Все неизвест- ные, кроме функции тока в конвективных членах уравнения переноса вихря (20), взяты c верхнего временного слоя (σ
6
= 1). Портрет мат- рицы соответствующей системы сеточных уравнений приведен на ри- сунке 4а. Ненулевые элементы отмечены звездочками (
∗
). Для сетки с числом узлов N
x
N
y
размерность матрицы линейной системы уравнений
2N
x
N
y
×
2N
x
N
y
•
Метод II. Уравнения (20) – (22) решаются совместно. В отличие от
Метода I, в этом подходе σ
6
= 0, т.е. на границе области вихрь, опре- деляемый по формулам Тома (22), вычисляется по значениям функции тока с предыдущего временного слоя. Портрет матрицы соответствую- щей системы линейных уравнений приведен на рисунке 4б. Ненулевых элементов в ней на 4(N
x
+ N
y
−
2) меньше, чем в аналогичной матрице для Метода I.

12
•
Метод III. Последовательный алгоритм: сначала из уравнений (20), (22)
при σ
6
= 0 определялся вихрь, затем по найденным значениям ω из уравнения (21) с граничными условиями ψ
∂Ω
= 0 вычислялась функ- ция тока. При таком подходе на каждом шаге по времени приходится решать две системы линейных уравнений размером N
x
N
y
×
N
x
N
y
. Порт- реты матриц приведены на рисунке 4. Заметим, что разностные схемы в алгоритмах II и III одинаковые, и эти два метода отличаются только способом обращения матрицы системы сеточных уравнений.
Описанные выше разностные алгоритмы реализованы в виде комплекса программ. Для решения систем уравнений использовалась библиотека от компании Intel: Intel Math Kernel Library, которая содержит эффективный параллельный алгоритм LU разложения матриц для решения систем ли- нейных уравнений с разреженными матрицами. Для хранения разреженных матриц использовался формат CSR (Compressed Sparse Row). В последова- тельном Методе III система уравнений (21) имеет матрицу с постоянными коэффициентами. Поэтому для решения систем уравнений (21) можно один раз выполнить LU факторизацию матрицы и в дальнейшем выполнять только обратный ход. При этом стоит отметить, что портреты матриц для уравнений
(20) и (21) одинаковы (рис. 4), то есть матрицу перестановок для них можно строить один раз.
0 5
10 15 20 25 30 35 40 45 50 0
5 10 15 20 25 30 35 40 45 50
а)
0 5
10 15 20 25 30 35 40 45 50 0
5 10 15 20 25 30 35 40 45 50
б)
Рис. 3. Портрет матрицы системы (20) − (22): а) для Метода I, б) для Метода II. N
x
=
N
y
= 4

13 0
5 10 15 20 25 0
5 10 15 20 25
а) Матрица для уравнения (20)
0 5
10 15 20 25 0
5 10 15 20 25
б) Матрица для уравнения (21)
Рис. 4. Портреты матриц для Метода III. N
x
= N
y
= 4.
5. Результаты расчетов
Описанные выше алгоритмы использовались для решения задачи об уста- новлении стационарного течения в конвекции Рэлея-Бенара при малой над- критичности. Рассматривалась прямоугольная область длиной L = 5, высоты
H = 1. На нижней границе поддерживается постоянная температура T
bot
= 1,
на верхней — T
top
= 0. В начальный момент времени (t=0):
V = 0,
T = 1
−
y + T
0
(x, y).
На линейный профиль температуры накладывалось возмущение
T
0
(x, y) = a sin (k
b
x) ,
a = 0.001 при y = H /3,
T
0
(x, y) = 0,
при y
6
= H /3.
Возмущение поля температуры порождает в слое жидкости движение в фор- ме валов с волновым числом k
b
. Рассматривалась малая надкритичность
Ra = 1800 при Pr = 1, то есть значение критического числа Рэлея превы- шено на 0.2%
−
0.6%Ra
cr
,
1
кинетическая энергия установившегося течения достигала E
kin
= 0.3.
1
Описание процедуры определения критического числа Рэлея приведено в приложении на странице 23.

14
Особенность подобных расчетов заключается в том, что возникающее те- чение развивается медленно. Длительность расчета может быть велика. Дви- жение будем считать установившимся, если кинетическая энергия E
kin
в те- чение 100 шагов по времени практически не изменяется,
|
E
j+100
kin
−
E
j
kin
|
<10
−10
,
где j — номер шага по времени.
5.1. Влияние аппроксимации граничных условий на величину
шага по времени
Исследовалась устойчивость рассматриваемых алгоритмов. Совместный метод I практически безусловно устойчивый [1], поэтому в расчетах шаг по времени был τ = 1t
ν
. Для определения условий устойчивости алгоритмов
II и III проведена серия расчетов для определения максимально возможного шага по времени τ
max
, при котором неустойчивость не возникает. Результа- ты представлены в таблице 1. Найденные значения τ
max для алгоритмов II
и III согласуются с известной теоретической оценкой, полученной методом
Неймана в случае одномерной линеаризованной модельной задачи [6]:
τ
max
∼
3 2
min
{
h
2
x
, h
2
y
}
(24)
Разностные схемы в методах II и III одинаковые, и отличаются способом решения соответствующих сеточных уравнений, что никак не сказалось на ограничении на шаг по времени.
Таблица 1. Шаг по времени τ
N
x
× N
y
Метод I
Метод II
Метод III
125 × 25 1
24.1E-4 24.1E-4 250 × 50 1
6E-4 6E-4 375 × 75 1
2.66E-4 2.66E-4
Основной причиной полученных ограничений на шаг по времени для Ме-
тодов II и III является то, что при использовании в граничных условиях (22)
вихря и функции тока с разных временных слоев нарушается дискретный аналог закона сохранения завихренности (9). Первоначально неустойчивость наблюдается в поведении вихря на границе области, а затем распространяет- ся внутрь.

15
5.2. Длительность расчета задачи
Проведена оценка затрат машинного времени для решения задачи до мо- мента выхода на стационар. В работе применялся параллельный LU алго- ритм, использующий перестановки ненулевых элементов в исходной матрице для минимизации заполненности матриц L и U. Процедуру перестановок можно выполнить один раз на первом шаге по времени, а затем использо- вать ее результат на всех последующих шагах. В таблице 2 показана дли- тельность всего расчета выхода на стационар на различных сетках. Время,
указанное в таблице 2, включает в себя время решения уравнений Навье-
Стокса и уравнения теплопроводности. При расчете использовались значения
τ
max найденные ранее. Таким образом, использование совместного алгоритма
Таблица 2. Длительность расчета [секунды]
N
x
× N
y
Метод I
Метод II
Метод III
125 × 25 11,21 344,80 243,79 250 × 50 73,69 7609,61 5828,0 375 × 75 181,51 54804,7 45137,8
с выбором функции тока в условии Тома с верхнего временного слоя показа- ло свою эффективность. Ограничения, накладываемые на шаги по времени в алгоритмах II и III, оказали существенное влияние на длительность расчета.
Ускорения, получаемые при использовании I алгоритма относительно II, ва- рьируются от 31 до 302 раз, а при использовании I относительно III получаем ускорения от 22 до 249 раз.
5.3. Затраты вычислительных ресурсов для расчета одного
шага по времени
Затраты машинного времени и затраты памяти на хранение матриц L и
U, необходимые для решения одного временного слоя, приведены в таблице
3. Время бралось средним по результатам 100 измерений. По этим парамет- рам совместные Метод I и Метод II одинаковы, так как в этих алгоритмах приходится обращать матрицу размера 2N
x
N
y
×
2N
x
N
y
. При этом обратить две матрицы размера N
x
N
y
×
N
x
N
y
в Методе III удается почти в два раза быстрее. Совместный алгоритм I дороже, чем Метод III: по времязатратам

16
почти в 2 раза, а по памяти — почти в 1, 3 раза. При этом для расчета всей задачи, за счет меньшего количества временных слоев, Метод I эффектив- нее, например, расчет одного конвективного времени t
ν
Методом I пройдет практически в 200 раз быстрее, чем Методом III.
Таблица 3. Затраты машинного времени и памяти для реализации одного шага по времени
Метод I
Метод II
Метод III
N
x
× N
y
время память время память время память
[секунды]
[Мбайт]
[секунды]
[Мбайт]
[секунды]
[Мбайт]
125 × 25 0,01 4,7 0,009 4,7 0,006 4,0 250 × 50 0,04 21,3 0,04 20,8 0,02 14,2 375 × 75 0,11 49,3 0,12 51,6 0,06 34,5
5.4. Зависимость времени выхода на стационар от величины
шага по времени. Расчеты Методом I
В таблице 4 показана зависимость длительности всего расчета и необхо- димое количество временных слоев для получения стационарного течения от шага по времени τ. С увеличением τ время расчета убывает, при τ > 16 вы- ход на стационар происходит практически за одно и тоже число временных слоев. То есть при τ > 16 τ приобретает смысл итерационного параметра,
а не шага по времени. Если необходимо исследовать не структуру стацио- нарного течения, а сценарий выхода на стационар, то расчет нужно вести с шагом τ < 1.
Таблица 4. Зависимость времени расчета от τ на сетке 375 × 75
Шаг
Время выхода
Количество
Шаг
Время выхода
Количество
τ
на стационар временных слоев
τ
на стационар временных слоев
0,1 296,94 2125 4
196,32 1261 0,5 240,43 1418 8
189,52 1250 1
225,41 1328 16 185,26 1244 2
200,70 1284 32 183,95 1242

17
5.5. Результаты распараллеливания. Затраты на один слой
по времени. Метод I
Все вычисления проводились с использованием четырехъядерного процес- сора Intel Core I7. В таблице 5 приведено время расчета одного шага по времени по Методу I на разном количестве ядер. Использование четырех ядер позволяет получить ускорение почти в три раза по сравнению с нерас- параллеленным LU методом.
Таблица 5. Длительность расчета одного шага по времени [секунды].
Количество
N
x
× N
y
ядер
125 × 25 250 × 50 375 × 75 1
0,0284 0,1336 0,3289 2
0,0131 0,0652 0,1814 4
0,0111 0,0482 0,1206
5.6. Влияние ширины ленты матрицы при совместном решении
уравнений Навье–Стокса
Использование совместного алгоритма решения уравнения Навье–Стокса связано с необходимостью решать систему уравнений с матрицей ленточной структуры (рис. 4а), расстояние между главной диагональю и самой дальней побочной диагональю равно 2N
x
или 2N
y
, в зависимости от порядка упорядо- ченности (ω
ij
, ψ
ij
) в векторе неизвестных. В работе проводилось исследование влияния ширины ленты матрицы на время расчета одного слоя по времени
Методом I. Рассматривались сетки с общим числом узлов N
x
N
y
= 160000,
распределенных по направлениям следующим образом: 400
×
400, 200
×
800,
800
×
200, 100
×
1600 и т.д.. Вектор неизвестных упорядочивался так, что ширина ленты равна 2N
x
. Результаты исследования приведены в таблице 6.
Полученный результат легко предсказуем: чем шире лента в исходной мат- рице, тем дольше идет расчет. В случае N
x
>> N
y
процедура перестановок ненулевых элементов в исходной матрице для минимизации заполненности матриц L и U , предусмотренная в библиотеке MKL, приводит к матрице с более узкой лентой. Поэтому длительности расчетов на сетках N
x
×
N
y
и
N
y
×
N
x
практически одинаковы.

18
Таблица 6. Длительность расчета одного шага по времени Методом I [секунды].
Сетка 400 × 400
Время 1, 30
N
x
× N
y
время
N
x
× N
y
время
200 × 800 1,25 800 × 200 1,25 100 × 1600 1,04 1600 × 100 1,00 50 × 3200 0,76 3200 × 50 0,65 25 × 6400 0,64 6400 × 25 0,52
Использование процедуры перестановок из библиотеки MKL позволило вести расчеты на подробных сетках, которые приводят к решению линейных систем сеточных уравнений большой размерности. Без предварительного по- строения матрицы перестановок не удается разместить матрицы L и U в оперативной памяти персонального компьютера. Например, для рассматрива- емой задачи с такой проблемой столкнулись на сетке 75
×
75. В таблице 7
приведены затраты машинного времени и требования по памяти для реше- нии системы уравнений размерности (2
·
50
·
50)
×
(2
·
50
·
50), возникающей в
Методе I на сетке 50
×
50. Если использовать алгоритм с предварительными перестановками, то расчеты на сетках 5
×
500 и 500
×
5 требуют практиче- ски одинаковое количество ресурсов. При этом расчет на этих сетках идет в два раза быстрее, чем на сетке 50
×
50, так как в матрицах, соответствую- щих сеткам 5
×
500 и 500
×
5, ширина ленты после перестановок уже. Без использования процедуры перестановок длительность расчета одного шага по времени определяется шириной ленты в разреженной матрице, именно поэтому самый долгий расчет — на сетке 500
×
5.
Таблица 7. Затраты машинного времени и памяти для реализации одного шага по времени.
Без предварительных перестановок
C предварительными перестановками время память время память
N
x
× N
y
[секунды]
[Мбайт]
[секунды]
[Мбайт]
50 × 50 0,04634 10,5 0,00924 4,4 5 × 500 0,00805 2,0 0,00452 2,1 500 × 5 1,86652 66,6 0,00390 2,0
Созданный комплекс программ позволяет эффективно вести расчеты на подробных сетках. Благодаря этому удалось детально разрешить структуру

19
течения для конвекции Рэлея–Бенара в прямоугольной области с жесткими боковыми границами (рис. 5). Обнаружить вторичные структуры течения в углах области (рис. 5б) [14].
-0,7
-0,6
-0,5
-0,4
-0,3
-0,2
-0,1 0,2 0,4 0,6 0,8 1
0 0,2 0,4 0,6 0,8 1
а.)
1 E - 6 3 E - 6 5 E -
6 7 E - 6 9 E - 6 0,025 0,05 0
0,025 0,05
б.)
Рис. 5. Конвекция Рэлея–Бенара в прямоугольной области. а.) Функция тока. б.) Вторич- ные структуры в левом нижнем углу. N
x
= N
y
= 450.

20
6. Заключение
Проведен анализ эффективности использования прямых методов для ре- шения разностных уравнений Навье–Стокса в переменных "вихрь – функция тока". Работа алгоритмов иллюстрировалась на примере задачи о конвекции
Рэлея–Бенара в прямоугольной области в условиях малой надкритичности. В
созданном программном комплексе для численного исследования конвекции
Рэлея-Бенара решение линейных систем сеточных уравнений осуществля- лось с помощью библиотеки Intel Math Kernel Library, которая содержит параллельный алгоритм LU разложения матриц для решения систем линей- ных уравнений с разреженными матрицами. Процедура предварительных пе- рестановок, реализованная в библиотеке MKL, позволила вести расчеты на подробных сетках. Показано, насколько совместный метод решения сеточ- ных уравнений Навье-Стокса дороже по затратам вычислительных ресурсов,
но за счет меньшего количества временных слоев, необходимых для решения рассматриваемого класса задач, использование совместного алгоритма дает ускорение расчета более чем в 200 раз по сравнению с последовательными методами решения уравнений Навье-Стокса.

21
Список литературы
1.
Мажорова О.С. Разностные методы решения уравнений динамики вяз- кой несжимаемой жидкости. Энциклопедия низкотемпературной плазмы.
Том VII-1. Математическое моделирование в низкотемпературной плазме.
Часть 2. Под ред. Попова Ю.П. М.: Янус-К. 2008. С. 235–253.
2.
Роуч П. Вычислительная гидродинамика.
Москва: Эдиториал УРСС,
1999. 247 c.
3.
Мажорова О.С., Попов Ю.П. Об одном алгоритме численного реше- ния двумерных уравнений Навье – Стокса // Препринты ИПМ им.
М.В.Келдыша. 1979. № 137. 24 c.
4.
Люмкис Е.Д. Об увеличении шага во времени при интегрировании урав- нений Навье – Стокса в переменных вихрь – функция тока. // Диффе- ренциальные уравнения. 1985. Т. 21, № 7. С. 1205–1217.
5.
Мажорова О.С., Попов Ю.П. О методах численного решения уравнений
Навье – Стокса. // Ж. Вычисл. Матем и мат. физики. 1980. Vol. 20,
no. 4. Pp. 1005 – 1020.
6.
Ермаков С.В., Мажорова О.С., Попов Ю.П. Математическое моделиро- вание задач электрофоретического разделения биосмесей. Ч.II. // Диф- ференц. уравн. 1992. Т. 28, № 12. С. 2129 – 2137.
7.
Флетчер К. Вычислительные методы в динамике жидкости.
Москва:
Мир, 1991. т.1. 62 с., т.2. 552 с.
8.
Ермаков М. К. Исследование возможностей матричных методов для ре- шения уравнений Навье – Стокса. // Физико-химическая кинетика в газовой динамике.
2010.
Т. 9.
http://chemphys.edu.ru/issues/2010–
9/articles/150/.
9.
Гершуни Г.З., Жуховицкий Е.М. Конвективная устойчивость несжимае- мой жидкости. Москва: Наука, 1972. 392 с.
10. Arakawa A. Computational design for long-term numerical integration of the equation of fluid motion: Two dimensional incompressible flow. // J.
Comput. Phys. 1966. Vol. 1. Pp. 119 – 143.

22 11. Моисеенко Б.Д., Фрязинов И.В. Консервативные разностные схемы для уравнений несжимаемой жидкости в переменных Эйлера. // Ж. Вычисл.
Матем и мат. физики. 1981. Т. 21. С. 1180 – 1191.
12. Thom A. The Flow Past Circular Cylinders at Low Speeds. // Proc.Roy.Soc.
London.Ser.A. 1933. Vol. 141. Pp. 651 –669.
13. Ковеня В.М., Яненко Н.Н. Метод расщепления в задачах газовой дина- мики. Отв. ред. Ю.И.Шокин. Новосибирск: Наука, Сибирское отделение,
1981. 304 с.
14. Ghia U., Ghia K.N., Shin C.T. High-RE solution for incompressible flow using the Navier-Stokes equations and multigrid method. // J. Comput.
Phys. 1982. Vol. 48. Pp. 387 – 411.

23
Приложение
В работе значения критических чисел Рэлея определялись численно. Для этого в статический профиль температуры вносилось возмущение вида
T
0
(x, y) = a sin (k
b
x) ,
a = 0.001 при y = H /3,
T
0
(x, y) = 0,
при y
6
= H /3,
где k
b
— волновое число возмущения, a — амплитуда. Если при рассматри- ваемом числе Рэлея в системе развивалось течение с кинетической энергией превосходящей 1.E
−
6, то считалось, что рассматриваемое число Рэлея лежит в области надкритичности. Значение числа Рэлея уменьшалось до тех пор,
пока все возмущения, задаваемые в статическом профиле, не начинали зату- хать, а жидкость возвращаться к статическому равновесию. За критическое число Рэлея принималось то значение, при котором энергия развивающегося течения менее 1.E
−
6, а его уменьшение на 1
−
2 приводит к затуханию любого возмущения. Таким образом были найдены следующие значения кри- тических чисел Рэлея для прямоугольной области длиной L = 5 и высотой H :
125
×
25 :
1770 < Ra
cr
< 1772 250
×
50 :
1775 < Ra
cr
< 1777 375
×
75 :
1776 < Ra
cr
< 1778

24
Содержание
1. Введение
3
2. Постановка задачи
5
3. Разностная схема
6
4. Аппроксимация по времени
11
5. Результаты расчетов
13
5.1. Влияние аппроксимации граничных условий на величину шага по времени . . . . . . . . . . . . . . . . . . . . . . . . . .
14 5.2. Длительность расчета задачи . . . . . . . . . . . . . . . . . . .
15 5.3. Затраты вычислительных ресурсов для расчета одного шага по времени . . . . . . . . . . . . . . . . . . . . . . . . . .
15 5.4. Зависимость времени выхода на стационар от величины шага по времени. Расчеты Методом I . . . . . . . . . . . . . .
16 5.5. Результаты распараллеливания. Затраты на один слой по времени. Метод I . . . . . . . . . . . . . . . . . . . . . . . .
17 5.6. Влияние ширины ленты матрицы при совместном решении уравнений Навье–Стокса
17
6. Заключение
20
Список литературы
21
Приложение
23