CC BY
136
УДК 517.958:536+519.62/64
Вестник СПбГУ. Сер. 10, 2013, вып. 4
Г. В. Кривовичев
ОБ ОДНОМ ВАРИАНТЕ МЕТОДА РЕШЕТОЧНЫХ УРАВНЕНИЙ БОЛЬЦМАНА
1. Введение. Метод решеточных уравнений Больцмана (lattice Boltzmann method, далее LBM) разработан для численного моделирования течений различных сред. Метод является альтернативой подходам, основанным на дискретизации уравнений, описывающих динамику сплошной среды [1—3]. При применении LBM решается система уравнений относительно функций распределения крупных частиц, из которой с помощью метода Чепмена-Энскога может быть получена система уравнений Навье-Стокса. Такие макроскопические величины как плотность, скорость, давление и температура приближенно вычисляются через значения функций распределения [1, 3].
Большая популярность LBM в настоящее время связана с простотой его алгоритма и широкими возможностями для распараллеливания вычислений. Метод практически идеально подходит для параллельной реализации на многопроцессорных графических ускорителях с использованием технологии CUDA [4-7], что является весьма актуальным [8, 9]. Кроме того, при реализации на многопроцессорных системах LBM демонстрирует увеличение производительности по сравнению с конечно-разностными методами решения уравнений гидродинамики и методом сглаженных частиц [10, 11]. Определенными преимуществами LBM обладает при расчетах течений со свободными поверхностями [11-13] и течений в пористых средах [14].
Данная работа посвящена описанию модифицированного варианта LBM, основанного на применении расщепления дифференциального оператора, входящего в уравнение Навье-Стокса, и идее мгновенной максвеллизации функций распределения, предложенной Б. Н. Четверушкиным [15, с. 46]. Эффективность модифицированного метода по сравнению с обычным LBM показана при решении известной тестовой задачи о течении в каверне (см. п. 4).
2. Метод решеточных уравнений Больцмана. В рамках LBM производится моделирование течения среды как динамики ансамбля крупных частиц с заданным конечным числом возможных скоростей. Область, в которой происходит течение, разбивается сеткой с ячейками определенной формы, что задает в ней так называемую решетку (lattice). За шаг по времени St частицы без взаимодействия друг с другом переходят между узлами решетки. Взаимодействие (абсолютно упругое соударение) может осуществляться только в узлах решетки [3].
В работе рассматриваются только плоские изотермические течения и решетка с ячейками в форме квадрата со стороной длины l. В этом случае можно использовать набор скоростей D2Q9 [1, 3]: Vj = Vvi, i = 1,..., 9, где V = l/St, а векторы Vj задаются следующим образом:
vi = (0,0), V2 = (1, 0), V3 = (0,1), V4 = (-1, 0), V5 = (0, -1), V6 = (1, 1), V7 = (-1, 1), V8 = (-1,-1), V9 = (1, -1).
Система разностных уравнений (так называемых решеточных уравнений Больц-мана (lattice Boltzmann equations, далее LBE)), описывающая динамику ансамбля
Кривовичев Герасим Владимирович — кандидат физико-математических наук, доцент, 199034, Санкт-Петербургский государственный университет; e-mail: gera1983k@bk.ru.
© Г. В. Кривовичев, 2013
крупных частиц, имеет вид [1-3]
/<(*,- + <Й,гы + У«Й) - = - /Р^г*,)) , (1)
где гк; = (жк, у;) - радиус-вектор узла решетки; ^ - узел временной равномерной сетки, построенной с шагом г = 19, — одночастичные функции распределения
крупных частиц со скоростями ; г - безразмерный параметр релаксации; -
функции, аппроксимирующие равновесные функции распределения Максвелла.
Система (1) может быть получена разными способами - как обобщение модели решеточного газа [3] или посредством дискретизации кинетического уравнения Бхатна-гара-Гросса-Крука [16].
В этой статье ограничимся рассмотрением течений вязкой несжимаемой жидкости. Тогда при проведении численных расчетов нужно находить только такие макроскопические характеристики как скорость среды и и давление р. При применении LBM эти величины вычисляются по формулам [1, 3]
9 9
р(г3, гы) = , гы), р^з, гы )u(tj, гы) = У/ (гз, г к;), (2)
3, г к;), р(Тз, г к;)
г=1 г=1
\V2.
в которых с = ^ I
Функции / (еч) представляются в таком виде [1, 3]:
/ , о(У • и) , 9 (У • и)2 3 и2\
+ + (3)
где
4 9 '
г = 1,
г = 2, 3,4, 5, ¿ = 6,7,8,9.
Выражения вида (3) справедливы лишь для случая малых значений безразмерного параметра М = |и|/У, который в рамках LBM имеет смысл числа Маха [3]. Необходимо отметить, что LBM, вообще говоря, предназначается для расчета течений сжимаемой среды. Случай малых М соответствует модели слабосжимаемой среды. Такая модель может с успехом использоваться для расчета течений несжимаемой жидкости [15, с. 126; 17, с. 136]. Многими авторами отмечалось [2, 18-20], что LBM связан с таким известным в вычислительной гидродинамике методом расчета течений вязкой несжимаемой жидкости как метод искусственной сжимаемости [17, с. 136-140; 20, 21].
Как показано в [2, с. 940-944], из системы (1) с учетом (2) и (3) при применении метода Чепмена-Энскога может быть получена следующая система уравнений относительно функций р и и:
V • и + 0(М2) = 0, (4)
о 1
-^ + (Уи)-и=—Ур + г/Ди + 0(М3), (5)
от р
где р - плотность жидкости. Система (4), (5) с точностью до членов одного порядка малости с М2 и М3, в которых учитывается эффект сжимаемости, соответствует системе уравнений вязкой несжимаемой жидкости. Это обеспечивает при малых М близость
решений задач для (1) с вычислением р и и по (2) с решениями задач для уравнений гидродинамики. В дальнейшем будем иметь ввиду только случай малых значений М и члены соответствующего порядка малости в уравнениях (4), (5) будем опускать.
Параметр V, фигурирующий в (5), имеет смысл кинематического коэффициента вязкости, и он связан с т соотношением [2, с. 943]
(«>
Один из недостатков LBM обусловлен тем, что система (1) является условно устойчивой по начальным условиям [22-24]. На практике эта проблема может быть решена за счет использования специальных граничных условий для функций распределения [25], схем с несколькими параметрами релаксации [26] или неявных вариантов LBE [27]. К сожалению, два последних подхода приводят к увеличению времени проведения расчетов. При этом при применении неявных LBE возникают определенные сложности при распараллеливании вычислений.
В настоящей статье предложен вариант LBM, использующий явную разностную схему (1), который позволяет улучшать устойчивость метода.
3. Модифицированный вариант 1БЫ. Предлагаемый подход основан на введении в векторное уравнение (5) параметра VI, имеющего размерность кинематического коэффициента вязкости. Это делается просто посредством добавления к правой части (5) слагаемого VlДu и его вычитания:
д и 1
+ (Уи) • и= (>-гл)Ди + г/1Ди. (7)
Вводя обозначения А\(и,р) = — (Уи)-и— ^Ур + г/хДи и ^(и) = (г/ —г/ 1)Ди, перепишем систему уравнений, составленную из (4) и (7):
ди
— =А1(и,р) + А2(и), V • и = 0. (8)
Следуя методу расщепления Н. Н. Яненко [17], предлагается решать задачи для (8) на промежутке времени [tj,tj+l], где tj+l = + 5Ь, таким образом: вводить промежуточный временной слой tj+l/2 = tj + 5Ь/2 и на промежутке [tj,tj+l/2] решать задачу для системы
ди
— =А1(и,р), V • и = 0, (9)
а затем, предполагая, что
p(tj+1, Гы) « p(tj+l/2, Гы), (10)
решать на промежутке [^+1/2,^+1 ] задачу для системы
(11)
для которой в качестве начального условия при Ь = tj+l/2 используется значение и, найденное при решении задачи для (9). При расчете на следующем временном промежутке как начальное условие для (9) применяется и, определенное при решении (11), а величина р задается в соответствии с (10).
Как можно видеть, предположение (10) позволяет избежать необходимости в расчете давления на промежутке [^3+1/2,^3+1 ]. В принципе использования (10) можно избежать, если по найденному при решении (11) полю скоростей в момент tj+l производить вычисление давления, решая задачу для уравнения Пуассона (например, см. [28, с. 276-278]). К сожалению, последнее привело бы к дополнительным ресурсоемким вычислениям. Предположение (10) является ни чем иным как упрощением, позволяющим избежать этой процедуры. Естественно, оно будет накладывать определенные ограничения на пределы применимости метода. Заметим, что в представленном в данной работе примере предположение (10) не сказывается на точности получаемых результатов.
При проведении численных расчетов вместо задачи для (9) предлагается решать задачу для системы (1) с параметром релаксации 71, вычисляемым через VI с использованием (6). Как известно [22, 23], посредством выбора параметра релаксации можно влиять на устойчивость схемы (1). В связи с этим за счет выбора величины VI, через которое рассчитывается 71, можно осуществлять стабилизацию системы (1).
На значения параметра VI накладывается условие VI > 0 (поскольку VI фигурирует в векторном уравнении движения на месте кинематического коэффициента вязкости, который неотрицателен), что влечет: 71 > 1/2. Потребуем также выполнения условия V—VI > 0, так как система (11) представляет собой систему из двух линейных уравнений диффузии, в которых V — VI имеет смысл коэффициента диффузии. Положительность этого коэффициента гарантирует устойчивость решений (11) по начальным условиям.
В настоящей работе предлагается выбирать следующее значение параметра релаксации 71: 71 = 1. Такой выбор обусловлен тем, что для стационарных течений оно соответствует внутренним точкам областей устойчивости в пространстве параметров [22-24, 29]. Единичное значение 71 соответствует и идее мгновенной максвеллизации функций распределения, предложенной Б. Н. Четверушкиным [15, с. 46], успешно применяемой при реализации кинетически-согласованных разностных схем решения задач газовой динамики. Последний момент поясним более подробно - перепишем (1) на промежутке \Ьз+ 5Ь/2\ при 71 = 1:
+ (12)
Система (12) описывает перенос значений функций вдоль характеристик системы уравнений Бхатнагара-Гросса-Крука (она представлена, например, в работе [16] на с. 242), выходящих из точек (Ьз, г к;). При использовании (12) предполагается, что при Ь = Ьз величины / в узлах (Ьз, г к;) совпадают со значениями / (е<. Таким образом, при проведении расчетов полагается, что в момент Ь = Ьз в узлах (3, г к;) происходит мгновенная максвеллизация функций распределения, поскольку, как указывалось выше, /(еч) аппроксимируют функции распределения Максвелла.
Из того, что 71 = 1, вытекает ограничение на параметр 7. Так как V — Vl > 0, и, согласно (6), Vl = I2/65Ь, получаем
¡2
Последнее подразумевает, что предлагаемый вариант LBM не будет подходить к задачам моделирования течений, для которых 7 « 1/2, что характерно, в частности, для турбулентных течений [30].
Система (11) представляет собой систему из двух линейных уравнений диффузии для компонент вектора и. Эти уравнения никак не связаны друг с другом, и задачи для них могут решаться одновременно и независимо друг от друга, что особенно важно для параллельной реализации расчетов на многопроцессорных системах. Для решения данных задач можно воспользоваться любой из разработанных и представленных в литературе (например, см. [31, 32]) разностных схем. При использовании различных схем можно добиться существенного улучшения устойчивости метода и смягчения условий на значения шага по времени.
Таким образом, применяя модифицированный вариант, можно улучшать устойчивость LBM. Это делается, во-первых, за счет выбора единичного значения Т1, которое принадлежит области устойчивости LBM. Во-вторых, устойчивость можно улучшить благодаря выбору разностной схемы решения задач для системы (11).
При программной реализации предложенного метода, чтобы не нарушать общей логики программы и не реализовывать принципиально другие алгоритмы, для решения уравнений из (11) можно воспользоваться одним из вариантов LBM для решения уравнения диффузии. Соответствующие явные схемы были предложены, например, в работах [33-35]. При этом, поскольку используются однотипные разностные системы на однотипных шаблонах, параллельная реализация алгоритма предлагаемого варианта LBM не будет отличаться от случая алгоритма обычного LBM.
В данной работе задачи для уравнений диффузии численно решались с помощью варианта LBM, предложенного D. А. Wolf-Gladrow [35], в котором используется малое число возможных скоростей из набора D2Q4: = Vу^, г = 1,..., 4, где
VI = (1, 0), У2 = (0,1), Уз = (-1, 0), У4 = (0,-1). Уравнение диффузии в узлах решетки решается следующим образом:
4
c(tj , Гк1) = ^ 9г (tj , Гк1),
¿=1
где g¿ - функции распределения, эволюция которых описывается системой, аналогичной (1), в которой равновесные распределения д(е9) физического смысла не имеют и вычисляются по формуле
Связь коэффициента диффузии Б с параметром т дается выражением [35]
V 2) 2бг
Для решения уравнений из (11), помимо явной схемы [35], применялась и полунеявная схема, основанная на методе переменных направлений [36]. Алгоритм этой схемы допускает естественное распараллеливание - при переходе с одного слоя на другой можно задавать столько вычислительных процессов, сколько узлов содержится в сетках по каждой из пространственных переменных [36, с. 465-470]. При решении тестовой задачи будет показано, что использование такой схемы позволяет существенно смягчить ограничение на шаг по времени, обеспечивающее устойчивость.
4. Решение задачи о течении в каверне. Применим предложенный вариант LBM к решению известной тестовой задачи вычислительной гидродинамики о плоском
течении в квадратной каверне [37, 38] и сравним результаты расчетов с полученными с помощью обычного LBM.
В задаче о течении в каверне рассматривается область П, границы которой параллельны осям декартовой прямоугольной системы координат:
П = {(x,y)\x е [0,P],y е [0,P],P> 0}. На границах ставятся только кинематические условия - условия прилипания [38]
ux(t, x, 0) = uy(t, x, 0) = 0, ux(t, x, P) = Uo, Uy(t, x,P) = 0, x е [0, P],
ux(t, 0,y) = uy(t, 0,y) = ux(t,P,y) = uy (t,P,y) = 0 У е [0,P), здесь U0 = const = 0.
В случае модифицированного LBM в качестве граничных условий для системы (12) выступали равновесные условия - функции распределения на границах приравнивались значениям f (eq). Для расчетов с использованием обычного LBM для системы (1) применялись условия из работы [39], для которых, как показано в [25], LBM имеет повышенную устойчивость по сравнению с другими известными типами граничных условий. При решении задач для системы (11) с помощью LBM использовался подход к реализации условий Дирихле, предложенный в [35]. Расчеты для двух вариантов LBM производились при помощи оригинальной программы, написанной на языке C+—+ (применялась некоммерческая среда разработки программного обеспечения Dev-C+—+ 4.9.9.2).
При проведении расчетов задавались значения такого критерия подобия как число Рейнольдса, в случае ньютоновской жидкости вычисляемого как
(13)
где L - линейный размер области, в которой происходит течение; U - характерная скорость. При решении задачи о течении в каверне в качестве L естественно рассматривать длину границы P, а U - величину Uo [38].
В настоящей работе автор не ставил целью проведение расчетов течений реальных сред, а просто решает тестовую задачу. Поэтому параметр v, через который вычисляется т, не задавался, а рассчитывался с использованием формулы (13) по заданным Re, U и L. Рассматривались малые и умеренные значения Re: 0.01, 1, 50, 100 при следующих параметрах: P = 1 м, U0 = 0.01 м/c. Изучался промежуток времени от 0 до 1200 с. Расчеты производились на сетках с одинаковым шагом по обеим декартовым координатам, с числом узлов, равным 50 х 50, 100 х 100, 150 х 150 и 200 х 200 соответственно.
В ходе расчетов находилось приближенное наибольшее значение шага по времени St, при котором не возникали численные неустойчивости для случая как обычного LBM, так и его модифицированного варианта. Ниже такое значение St обозначим через St.
На рисунке для сравнения представлены графики численных решений, отнесенных к Uo, полученных с помощью модифицированного LBM в случае сетки из 200x200 узлов и данных, взятых из работ [37, 38, 40]. На рисунке, а приведены графики ux/Uo в точках прямой {x = 0.5P, y е [0, P]}, а на б - графики uy/Uo в точках прямой {x е [0, P], y = 0.5P}. Как можно видеть, результаты, полученные с помощью модифицированного LBM, хорошо согласуются с литературными данными.
иу/и0
-0.2 0 0.2 0.4 0.6 0.8 1
их/и0
-0.4-0.2 0 0.2 0.4 0.6 0.8 1
0.25
иг/ип
III
-0.4-0.2 0 0.2 0.4 0.6 0.8 1
их/Щ
"у
0.2
0.1
0 -0.1
-0.2 -0.3
и„/и0
о
1 х
0.2 0.4 0.6 0.8 1 х
0.2 0.4 0.6 0.8 1 х
Графики компонент вектора скорости при Ев = 1 (I), 50 (II) и 100 (III) 1 — результаты расчетов по модифицированному ЬБМ; 2 — результаты из [40] (I), [37] (II), [38] (III).
Объяснение а и б в тексте.
Для исследования влияния разбиения сетки на точность результатов производилось вычисление следующей величины:
1=^РХ+Ру (14)
В (14) 1Х и 1у имеют смысл среднеквадратичных отклонений полученных результатов от представленных в работах [37, 38, 40]:
1 М 1 М
Р = — - их(0.5Р,у^)2, 1у = — 53(г4у(Т,х<,0.5Р) - 0.5Р))2,
i=0 I=0
где Т - длина рассматриваемого временного промежутка; их и иу - взятые из литературы значения компонент вектора скорости в узлах равномерной сетки, в которой отрезок [0,Р] разбивается на N частей. При вычислении 1Х и 1у бралось N = 10.
В табл. 1 представлена величина I в случае использования схемы D. А. Wolf-Gladrow: как можно видеть, она уменьшается при измельчении сетки для случаев Ев = 1, 50 и 100, что свидетельствует о практической сходимости модифицированного LBM.
Таблица 1. Значения I для Кв = 1, 50,100 и различных разбиений сетки в случае модифицированного ЬБМ при использовании схемы из [35] решения задачи Дирихле для системы (11)
Значение Re 50 X 50 узлов 100 X 100 узлов 150 X 150 узлов 200 X 200 узлов
1 1.042 • Ю-4 2.263 • КГЬ 8.077 • КГЬ 4.298 • 10~ь
50 9.646 • КГ4 5.067- КГ4 3.088 • КГ4 2.472 • 10~4
100 2.130 • КГ3 1.215 • КГ3 1.021 • КГ3 8.462 • 10~4
В табл. 2 приведены значения для случаев обычного LBM и его модификации при использовании схемы из [35]. Как можно заметить, значения для модифицированного LBM в случаях всех сеток больше в 2-4 раза, чем для обычного LBM. Этот факт позволяет говорить о том, что даже при применении явной схемы для решения (11) требование на шаг, связанное с обеспечением устойчивости, для модифицированного LBM является более мягким, чем в случае обычного LBM. К тому же модифицированный метод можно считать экономичным - при расчетах с более крупным шагом снижается размерность хранимых в памяти массивов данных.
Таблица 2. Значения для различных Кв и разбиений сетки при применении разных методов
Значение Re | 50 X 50 узлов | 100 X 100 узлов | 150 X 150 узлов | 200 X 200 узлов
Обычный LBM
0.01 5.555 • КГЬ 1.667- КГЬ 1.001 • КГЬ 4.166 • 10~ь
1 6.250 • КГ3 2.777- 10~л 1.724 • 10~л 1.190 • 10~л
50 1.497- Ю-1 6.822 • 10-'2 4.363 • 10-'2 3.188 • Ю-'2
100 1.667- КГ1 7.419 • 1Q-'2 4.675 • 1Q-'2 3.390 • 1Q-'2
Модифицированный LBM
0.01 2.502 • КГ4 8.196 • КГЬ 5.319 • КГЬ 2.083 • 10~ь
1 2.272 • 10-'2 5.548 • Ю-3 3.448 • Ю-3 2.503 • Ю-3
50 5.268 • КГ1 1.925 • КГ1 1.190 • КГ1 8.001 • 10-'2
100 7.824 • КГ1 3.573 • КГ1 1.689 • КГ1 1.149 • 1Q-1
При решении задачи Дирихле для системы (11) полунеявным методом переменных направлений удалось получить существенно большие значения St (табл. 3), например при Re = 1 они оказались больше примерно в 6 раз, при Re = 0.01 - приблизительно в 20 раз. Таким образом, посредством использования специальной схемы условие на шаг удалось смягчить еще сильнее, что приводит к значительной экономии вычислительных ресурсов. Это преимущество будет особенно существенно в практических приложениях, связанных с расчетом пространственных течений в областях сложной формы, когда используются сетки с миллионами узлов.
Таблица 3. Значения St для случаев Re = 0.01, Re =1 и различных разбиений сетки в случае модифицированного LBM при применении метода переменных направлений решения задачи Дирихле для системы (11)
Значение Re 50 X 50 узлов 100 X 100 узлов 150 X 150 узлов 200 X 200 узлов
0.01 1.102 • 10-á 3.250 • Ю-4 1.933 • Ю-4 8.011 • 10~ь
1 3.701 • 10-'2 1.670 • 10-'2 9.901 • КГ3 7.110 • 10~3
Поскольку при решении задач для системы (11) можно применять схемы, алгоритмы которых обладают естественным параллелизмом, модифицированный LBM обладает тем же преимуществом, связанным с удобством реализации алгоритма на многопроцессорных системах, что и обычный LBM.
5. Заключение. Полученные результаты позволяют прийти к выводу, что модифицированный LBM для задачи о течении в каверне оказался более предпочтительным из-за возможности проведения расчетов с большим значением шага по времени, что говорит о его экономичности и улучшенной устойчивости по сравнению с обычным LBM.
Как уже отмечалось, метод может использоваться только при малых и умеренных значениях Re. Возможные пределы применимости метода, связанные с предположением (10), еще предстоит выяснить. Но по всей видимости, метод будет хорошо подходить для расчетов установившихся течений в пористых средах, для которых характерно малое значение Re и предположение (10) является оправданным, поскольку давление слабо меняется во времени [41]. Заметим, что такие течения происходят в костях млекопитающих и мягких тканях человеческого организма [41]. Предложенный метод может быть распространен на случай пространственных течений и течений со свободной поверхностью.
Литература
1. Chen S., Doolen G. D. Lattice Boltzmann method for fluid flows // Ann. Rev. of Fluid Mech. 1998. Vol. 30. P. 329-364.
2. He X., Luo L. S. Lattice Boltzmann model for the incompressible Navier-Stokes equation // J. of Stat. Phys. 1997. Vol. 88, N 3/4. P. 927-944.
3. Wolf-Gladrow D. A. Lattice-gas cellular automata and lattice Boltzmann models — an introduction. Berlin: Springer, 2005. 311 p.
4. Бикулов Д. А., Сенин Д. С., Демин Д. С., Дмитриев А. В., Грачев Н. Е. Реализация метода решеточных уравнений Больцмана для расчетов на GPU-кластере // Вычислительные методы и программирование. 2012. Т. 13, № 1. C. 221-228.
5. Грачев Н. Е., Дмитриев А. В., Сенин Д. С. Моделирование динамики газа при помощи решеточного метода Больцмана // Вычислительные методы и программирование. 2011. Т. 12, № 1. C. 227-231.
6. Rinaldi P. R., Dari E. A., Venere M. J., Clansse A. A lattice Boltzmann solver for 3D fluid simulation on GPU // Simul. Model. Pract. and Theor. 2012. Vol. 25. P. 163-171.
7. Xiong Q. G., Li B, Xu J., Fang X. J., Wang X. W., Wang L. M., He X. F., Ge W. Efficient parallel implementation of the lattice Boltzmann method on large clusters of graphic processing units // Comput. Sc. and Technol. 2012. Vol. 57, N 7. P. 707-715.
8. Юдин И. П., Перепелкин Е. Е. Использование параллельных вычислений на графических процессорах при исследовании пропускной способности канала транспортировки пучка ионов с учетом пространственного заряда // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2012. Вып. 3. С. 103-112.
9. Раба Н. О. Разработка и реализация алгоритма расчета коагуляции в модели облаков со смешанной фазой с использованием технологии CUDA // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2011. Вып. 4. С. 94-104.
10. Velivelli A. C., Bryden K. M. Parallel peformance and accuracy of lattice Boltzmann and traditional finite difference methods for solving the unsteady two-dimensional Burger's equation // Physica A. 2006. Vol. 362. P. 139-145.
11. Гугушвили И. В., Евстигнеев Н. М. Некоторые результаты для различных методов моделирования несжимаемой гидродинамики свободной поверхностью на графических процессорах // Учен. записки. Электрон. науч. журн. Курск. гос. ун-та. 2010. Вып. 4. С. 15-23.
12. Schreiber M., Neumann P., Zimmer S., Bungartza H. J. Free-surface lattice-Boltzmann simulation on many-core architectures // Procedia Comput. Sc. 2011. Vol. 4. P. 984-993.
13. Zhao Z., Huang P., Li Y., Li J. A lattice Boltzmann method for viscous free surface waves in two dimensions // Intern. J. for Numer. Meth. in Fluids. 2013. Vol. 71. P. 223-248.
14. Pan C., Luo L. S., Miller C. T. Lattice Boltzmann model for incompressible flows through porous media // Phys. Rev. E. 2002. Vol. 66. P. 036304-1-036304-9.
15. Четверушкин Б. Н. Кинетически-согласованные схемы в газовой динамике. М.: Изд-во Моск. ун-та, 1999. 232 с.
16. Abe T. Derivation of the lattice Boltzmann method by means of the discrete ordinate method for the Boltzmann equation // J. of Comput. Phys. 1997. Vol. 131, N 1. P. 241-246.
17. Яненко Н. Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967. 197 c.
18. He X., Doolen G. D, Clark T. Comparison of the lattice Boltzmann method and the artificial compressibility method for Navier-Stokes equations // J. of Comput. Phys. 2002. Vol. 179. P. 439-451.
19. Asinari P., Ohwada T., Chiavazzo E., Rienzo A. F. Link-wise artificial compressibility method // J. of Comput. Phys. 2012. Vol. 231. P. 5109-5143.
20. Ohwada T., Asinari P. Artificial compressibility method revisited: asymptotic numerical method for the incompressible Navier-Stokes equations // J. of Comput. Phys. 2010. Vol. 229. P. 1698-1723.
21. Chorin A. J. A numerical method for solving incompressible viscous flow problems // J. of Comput. Phys. 1967. Vol. 2. P. 12-26.
22. ¡Sterling J. D., Chen S. Stability analysis of lattice Boltzmann methods // J. of Comput. Phys. 1996. Vol. 123. P. 196-206.
23. Кривовичев Г. В. Об устойчивости решеточной кинетической схемы Больцмана для расчета плоских течений // Вычислительные методы и программирование. 2011. Т. 12, № 1. С. 194-204.
24. Кривовичев Г. В. Исследование устойчивости явных конечно-разностных решеточных кинетических схем Больцмана // Вычислительные методы и программирование. 2012. Т. 13, № 1. С. 332-340.
25. Семенов С. А., Кривовичев Г. В. Численное исследование подходов к реализации граничных условий в методе решеточных уравнений Больцмана // Процессы управления и устойчивость: Труды 43-й междунар. науч. конференции / под ред. А. С. Еремина, Н. В. Смирнова. СПб.: Издат. дом С.-Петерб. ун-та, 2012. C. 196-201.
26. d'Humieres D., Ginzburg I., Krafczyk M., Lallemand P., Luo L. S. Multiple-relaxation-time lattice Boltzmann models in three dimensions // Philos. Transact. of Royal Soc. of London A. 2002. Vol. 360. P. 437-451.
27. Кривовичев Г. В. О применении интегро-интерполяционного метода к построению одноша-говых решеточных кинетических схем Больцмана // Вычислительные методы и программирование. 2012. Т. 13, № 1. C. 25-34.
28. Роуч П. Вычислительная гидродинамика / пер. с англ. В. А. Гущина, В. Я. Митницкого; под ред. П. И. Чушкина. М.: Мир, 1972. 612 с. (Roache P. Computational fluid dynamics.)
29. Nourgaliev R. R., Dinh T. N., Theofanous T. G., Joseph D. The lattice Boltzmann equation method: theoretical interpretation, numerics and implications // Intern. J. of Multip. Flow. 2003. Vol. 29. P. 117-169.
30. Martinez D. O., Matthaeus W. H., Chen S., Montgomery D. C. Comparison of spectral method and lattice Boltzmann simulations of two-dimensional hydrodynamics // Phys. of Fluids. 1994. Vol. 6, N 3. P. 1285-1298.
31. Самарский А. А., Вабищевич П. Н. Вычислительная теплопередача. М.: Книжный дом «ЛИБ-РОКОМ», 2009. 784 с.
32. Самарский А. А., Вабищевич П. Н. Численные методы решения задач конвекции—диффузии. М.: Книжный дом «ЛИБРОКОМ», 2009. 248 с.
33. Blaak R., Sloot P. M. A. Lattice dependence of reaction-diffusion in lattice Boltzmann modeling // Comput. Phys. Commun. 2000. Vol. 129. P. 256-266.
34. Wang J., Wang M., Li Z. Lattice Poisson-Boltzmann simulations of electro-osmotic flows in microchannels // J. of Coll. and Interf. Sc. 2006. Vol. 296. P. 729-736.
35. Wolf-Gladrow D. A. A lattice Boltzmann equation for diffusion // J. of Stat. Phys. 1995. Vol. 79, N 5-6. P. 1023-1032.
36. Киреев В. И., Пантелеев А. В. Численные методы в примерах и задачах. М.: Высшая школа, 2004. 480 c.
37. Chen S., Tolke J., Krafczyk M. A new method for the numerical solution of vorticity - streamfunction formulations // Comput. Meth. in Appl. Mech. and Engineer. 2008. Vol. 198. P. 367-376.
38. Ghia U., Ghia K. N., Shin C. T. High-Re solutions for incompressible flow using the Navier-Stokes equations and multigrid method //J. Comp. Physics. 1982. Vol. 48. P. 387-411.
39. Le Coupanec E., Verschaeve J. C. G. A mass conserving boundary condition for the lattice Boltzmann method for tangentially moving walls // Math. and Comput. in Simulation. 2011. Vol. 81, N 12. P. 2632-2645.
40. Семин Л. Г., Шапеев В. П. Метод коллокаций и наименьших квадратов для уравнений Навье-Стокса // Вычислительные технологии. 1998. Т. 3, № 3. C. 72-84.
41. Hussein M. A., Becker T. Numerical modelling of shear and normal stress of micro-porous ceramics for stimulated in-vitro cultivation of bone cells // Microfluid. and Nanofluid. 2010. Vol. 8. P. 665-675.
Статья рекомендована к печати проф. Н. В. Егоровым. Статья поступила в редакцию 30 мая 2013 г.
