CC BY
43
решения) к смежной точке, соответствующей улучшенному значению целевой функции. Максимальное количество итераций, необходимых для получения оптимума (т. е. количество базисных решений), не превосходит и!/[(и - т)!т!], где п - число переменных, а т - число уравнений модели линейного программирования, представленной в стандартной форме.
Предложенные методики определения оптимального парка с различными типами воздушных судов (с учётом условий их применения), маршрутов воздушного сообщения, позволяют рекомендовать авиапредприятиям применять воздушные суда с наибольшей эффективностью.
Работа выполнена в рамках Государственного контракта № П295 от 24.07.2009 ФЦП "Научные и научно-педагогические кадры инновационной России"
СПИСОК ЛИТЕРАТУРЫ
1. Бузулук В.И. Оптимизация траекторий движения аэрокосмических летательных аппаратов. М.: ООО "Марийское рекламно-издательское полиграфическое предприятие", 2008. 476 с.
2. Трусов П.В. Введение в математическое моделирование: Учеб. пособие. М.: Логос, 2004. 136 с.
3. Припадчев А.Д. Математические модели, применяемые для пассажирских перевозок. Системы проектирования, моделирования, подготовки производства и управление проектами CAD/CAM/CAE/PDM // Сб. ст. III Междунар. науч.-практ. конф. Пенза: Приволжский Дом знаний, 2009. С. 59-61.
4. Самарский А.А., Михайлов А.П. Математическое моделирование. Идеи. Методы. Примеры. М.: Физматлит, 2001. 320 с.
5. Седжвик Р. Фундаментальные алгоритмы. СПб.: ДиаСофтЮП, 2003. 672 с.
6. Арнольд В. И. Жесткие и мягкие математические модели. М.: МЦНМО, 2004. 144 с.
7. Альфред Ахо В., Хопкрофт Д., Джеффери Ульман Д. Структуры данных и алгоритмы. М.: Изд. дом "Вильямс", 2000. С. 384-390.
8. Безручко Б.П., Смирнов Д.А. Математическое моделирование и хаотические временные ряды. Саратов: ГосУНЦ "Колледж", 2005. 168 с.
9. Волков Е.А. Численные методы. М.: Физматлит,
2003. 298 с.
10. Исследование операций в экономике: Учеб. пособие для вузов / Под ред. Н.Ш. Кремера. М.: ЮНИТИ,
2004. 407 с.
УДК 519.63:004
Н.М. Евстигнеев
ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ УРАВНЕНИЙ НАВЬЕ-СТОКСА НА НЕСТРУКТУРИРОВАННЫХ СЕТКАХ С ПРИМЕНЕНИЕМ ЛАГРАНЖЕВО-ЭЙЛЕРОВОГО МЕТОДА
1. Исходные уравнения
Будем рассматривать систему уравнений Навье-Стокса, записанную в интегральном виде, и начально-краевую задачу для них как:
Пусть задана стационарная произвольная ограниченная область Ое Ш3 с границами ГОе Ш2. Найти вектор-функцию скорости ^О* [0,Г] ^ Ш3 и скалярную функцию давления Р:0х[0,7] ^ Ш такие, чтобы для любой замкнутой области внутри О они удовлетворяли интегральным соотношениям:
[дШ + ■ Ш - <¡1 -Ш-Р,- =0,(1) о эа за
где д = [0;и;у;^]г - вектор переменных количества движения, V = в декартовой системе координат; ■ п =[(©);(и© + пхР); (у© + п Р);
+ пг Р)]г - невязкие потоки для произвольной подобласти внутри О; п = {п, п, и } - вектор нормали к границам замкнутой подобласти;
?! ■ п = [0;(п т + пх + пт + т + т ); (п т + пт +
4 х д^ у ху г хг хи гх 7 \ х ух у уу
+ пт + т + т ); (пт + пт + пт + т + т )]г-
г уг ух у2/~\х гх у гу г гг гх гу' J
вектор вязких потоков для произвольном подобласти; 0 = nu + nv + nw; {ij} = {1, 2, 3} -1 ' ' Э Vj
декартовы направления; т, = —(—-*- + —- тен-
R axj axt
зор вязких напряжений; g - единичный вектор направления действия гравитации; R = UL/v, Fr = U/^L - числа Рейнольдса и Фруда, соответственно, составленные по характерному масштабу течения; m^ - тензор вязких напряжений, обусловленный наличием оператора произвольного осреднения для турбулентного течения (определение дано в п. 4); Fs - функция, отвечающая за поверхностное натяжение в случае наличия фазы раздела сред (определение дано в (27) п. 5). В данной статье ограничивается рассмотрение фазы раздела жидкость-атмосфера (течение со свободной поверхностью).
2. Метод решения
Существует достаточно обширное семейство различных методов дискретизации О и решения НКЗ для (1) выбранной дискретизации [1]. Обзор указанных методов выходит за рамки этой статьи. Остановимся на методе конечного объёма, на неструктурированной сетке, поскольку данный метод позволяет строго выполнять балансовые соотношения (1), является консервативным и позволяет применять методы повышенного порядка [10]. Применение неструктурированной сетки обусловлено рассмотрением О сложной произвольной топологии, встречающейся в большинстве практически важных задач.
Для согласования поля скалярной функции давления и соленоидальности по вектор-функции скорости используется метод расщепления, впервые предложенный академиком О.М. Белоцерков-ским, основанный на методе Маркеров и Ячеек (MAC) [1]. Применяя оператор проекции скалярной функции давления на оператор дивергенции и расщепляя (1) по физическим процессам, получим для метода конечного объёма [10]: Шаг 1. Невязкие потоки и внешние силы:
и+1/4 П
и, - и i
At
W =
К-1 3
= -Z5>V-Ы-Fs-Fr-'g. (2Л)
k=Qj=l
Ш а г 2. Вязкие потоки:
п+2/4 л+1/4 и, - u,
к 3 1 ?и
k=0 j=1 К dXi
At
W =
п+2/4
+ V2/4)/<t-ny>t.5t].(2.2)
Шаг 3. Уравнение Пуассона для согласования поля давления при условии что div(Fn+1) = 0 . V2P = - VVn+2/4/At. (2.3)
Ш а г 4. Корректор соленоидальности вектор-функции скорости: Vn+1 = Vn+2/4 - At VP. (2.4) Здесь: i, j = {1, 2, 3} - декартовые направления; n -шаг по времени; k = [0..^) - количество сторон элемента на неструктурированной сетке; At - шаг по времени; W - объём элемента сетки; Sk - площадь k-й стороны элемента; n - составляющая вектора единичной нормали грани k по декартовому направлению j. Оператор набла (V) также записывается методом конечного объёма.
Для (2) ставится следующая начально-краевая задача:
V(t, Г1П) = К(Г1П); д V(t, Г2П)/5(п/ Г2О) = = V0 = const; V(0, О) = V0,
div(V0) = 0; P(t, [ГОщ]) = др = = const; P(t, [ГОсг|]) = P0 = const , P(0, Q) = P0. (3)
Где п - произвольная гиперплоскость контакта жидкость-атмосфера.
Исходя из рассмотренного метода расщепления очевидно, что устойчивость всего метода определяется исходя из шагов 1 и 2, поскольку шаг 4 является безусловно устойчивым.
Рассмотрим отдельно каждый шаг и способ его решения.
3. Численный Лагранжево-Эйлеров метод для невязких потоков
Шаг 1 для схемы расщепления можно переписать в виде:
рП+1/4 _ рп. _
W Л
= "Т(4)
к=ом
где F - скалярная функция (аналог декартовой составляющей V (x,t)).
Разобьём, для примера, исходную область О 3-симплексами (тетраэдрами) общим количеством N, и будем решать задачу (3) для уравнения (4), переписанного в следующем виде [2]:
— = 0; —=V(x,t);xe Q[0. JV]с9?3. (5) dt dt
Алгоритм для решения (5) записывается следующим образом:
1. Лагранжев этап. Каждая точка в О, связан -ная с центром массы тетраэдра с координатами xi,
—*
переносится в координаты х на момент времени t (в предыдущее положение) в О путём интегрирования траектории:
<+д«
х*1=х,- |у(т)Л. (6)
г
2. Эйлеров этап. Нахождение элемента (локализация в О) в котором находится точка с координатами х . и вычисление значения скалярной величины Е путём интерполяции соседних величин Е из тетраэдров, являющихся соседями К-го порядка:
е( х >=1(Е{ % ¿а (7)
где I - оператор интерполяции; х К - координаты соседей вплоть до соседей порядка К.
3. Регуляризационный этап. Скалярное поле обновляется на шаге по времени t +Д^
Е( х ^ + Дt) = Е( х * (8)
Процедура (6)-(8) при условии выполнения условия ТУВ [9,10] является монотонной, консервативной (т. е. соответствующей записи (4)) [5] и:
Утверждение 1. Процедура решения (5) с помощью метода (6)-(8) - абсолютно устойчива по временному шагу.
Доказательство. Рассмотрение линейного приближения (4) для одномерного случая позволяет искать решение как суперпозицию Фурье гармоник вида:
Е = а-ехр[гк(х - V)], где к = 2я/(№ W/V). (9)
Переписывая: Е = а-ехр(-/юОехр(/кх), с учётом ю = (9) перепишется для дискретного разбиения (Дх = W/ V) как:
Е^ = аАп ехр(гктДх); Ап = ехр(-г'юД0. (10)
Координата х*, определяемая по (7), в одномерном случае находится как: х = х - VДt, и тогда -х + х * = - VДt, ^ СЕЬ = VДt/Дх, где СЕЬ - число Куранта-Фридрихса-Леви. Введём а=CFL - [СЕ£] и, переписав интерполяционную процедуру (7) для К = 1, получим решение (8) в виде:
Е(х ^ + ДО = аЕ(х ^ + (1 - а) Е(х^), (11)
где ае [0,1]<=Ш, р = [CЕL]cZ. Тогда, подставляя (10) в (11):
а-Ая+1 ехр(г'кЖх) = ааАп ехр(гк(т - p - 1)Дх) + +(1-а)аАяехр(/к(т -р)Дх), или: А = аexp(/k(-p - 1)Дх) + (1 - а) ехр(гк(-р)Дх), ^ А = ехр(-г кр Дх)[(1 - а) + а ехр(-г кДх)] (12)
Возведём (12) в квадрат по модулю:
\a\2= \exp(-ikpAx)\2^\[(1 - a) + a exp(-z'Mx)]|2 =
12
= \ (1 - a) + a cos kAx - i a sin kAx\ , ^
\a\2= [(1 - a) + a cos kAx]2 + [a sin kAx]2 = = (1 - a)2 + 2(1 - a)a cos kAx + a2 cos2 kAx + + a2sin2 kAx = (1 - 2a + a2)+2a(1 - a)cos kAx + a2 = = 1 - 2a(1 - a)[1 - cos kAx]. (13)
Поскольку в (13) 0 < (1- cos kAx) < 2, то максимальное значение в (13) имеет вид:
\a\2= 1 - 4a + 4a2 = (1 - 2a)2 < 1, (14)
а минимальное значение определяет \a \2 как:
\a\2= 1. (15)
Из (14) и (15) следует, что решения, описываемые (9) не будут расти при переходе с одного шага по времени на другой.
Таким образом, основными этапами решения задачи (5) являются нахождение траектории по выражению (6); определение номера элемента (3-симплекса), в котором находится точка с вычисленными по (6) координатами; интерполяция F то (7). Более детально данные этапы применительно к сетке тетраэдров рассмотрены в [18]. В данной статье кратко описаны основные особенности.
Процедура регуляризации элемента, достаточно просто выполняющаяся для структурированной дискретной сетки разбития исходной области [4], является значительной трудностью при переходе на неструктурированную сетку (например, тетраэдров), поскольку просмотр всех ячеек приводит к трудозатратам порядка N2, что для задач с большим количеством элементов недопустимо. Существует несколько ускоренных процедур локализации, т. е. определения точного дискретного элемента с координатами {x*; y*; z*.). Очень хороший обзор по данным процедурам для неструктурированных сеток дан в [4], где показано, что наиболее точным и быстрым методом является метод, предложенный в [3]. Детально не останавливаясь на всех методах, здесь используется метод, аналогичный [3], только более быстрый и устойчивый метод поиска и локализации. Процедура представляет собой алгоритм направленного поиска, который определяет грани симплекса, через которые частица выходит из него, для чего используется условие "частица слева от грани". Далее алгоритм проходит по всем элементам, где частица пересекла грань, пока частица
не окажется справа от грани для всех граней симплекса. Алгоритм достаточно устойчив и с вычислительной точки зрения имеет производительность, пропорциональную Л'С'М, где С -количество элементов, пересёкших траекторию, Л - количество граней симплексов (четыре для тетраэдров, пять для пирамид и т. д.).
Для выполнения интерполяции используется диаграмма Вороного на системе симплексов. Нахождение значения скалярной функции Е в точке х ведётся с применением интерполяции в ячейке политопа Вороного, вершинами которого являются центры масс прилегающих симплексов. Так, для тетраэдра t вес в значении Ех записывается как взвешенная сумма значений вершин политопа. Вес вершины t со значением Е вычисляется как:
МГ,(х) =
(16)
при этом:
т
2>(х),=1. (17)
г=\
Здесь о - множество граней политопа, примыкающих к вершине V, п\ - нормаль грани; х - расстояние от точки до грани; ё - смещение плоскости грани; т - количество граней в политопе ячейки Вороного, образованного системой тетраэдров; - детерминант матрицы единичных нормалей граней в ог Если сетку тетраэдров для произвольной области строить с использованием алгоритма Делоне [11], то выражение (16) значительно упрощается, поскольку выполняются два основных свойства: в сетке Делоне вектора, соединяющие вершины, всегда перпендикулярны граням диаграмм Вороного; объём тетраэдра t определяется как где Et - матрица, сфор мированная тремя векторами сторон тетраэдра е 1, е2, е3, имеющими общую точку из числа вершин тетраэдра t. Тогда выражение (16) может быть переписано как:
Щ(х) =
ш,
I е1 II е2
1ЛГ, I
(18)
6 УоЩ)
где: Vol(t) - объём тетраэдра V, е' - вектор грани из точки, принадлежащей политопу Вороного
к другим вершинам % с - центр тяжести % х - координаты точки, для которой производится интерполяция.
После вычисления ^(х) по (18) и поскольку справедливо (17), Ех определяется как:
(19)
которое используется на Регуляризационном этапе. На данном этапе учитываются и граничные условия (4) очевидным образом - для условия Дирихле вычисляется виртуальная ячейка со сходной топологией, что и соседи и ставится явное значение скалярной функции (4). Для условия Неймана значение грани при умножении в знаменателе заменяется значением производной скалярной функции по направлению.
При вычислении траектории в п.2 применяется интерполяция скорости на шагах по времени t и t+Дt. Но при этом, в случае немонотонности функции Е, могут возникнуть нефизические осцилляции. Для компенсации данного эффекта и для повышения порядка аппроксимации используется ТУБ алгоритм совместно с алгоритмом BFECC [12] (компенсация и коррекция ошибки прямым и обратным ходом). Более детально [6]. Получаемая схема имеет второй порядок точности по пространству и по времени.
4. Метод конечного объёма для вязких потоков.
Турбулентные напряжения, модель Динамики Больших Вихрей
Метод конечного объёма для вязких потоков записывается через перекрёстное интегрирование на гранях симплексов. Так, для простоты изложения опуская тензор т.. , получим: At к 3 7)и
1 (20)
л+2/4 л+1/4
и, - и.
Переписывая (20) с учётом топологии симплексов, получим:
л+2/4 л+1/4 и, - и, :
А; 1
л+2/4 л+2/4 и. - и,
п1,к Як
«м5*.(21)
где расстояние от центра хранения дискретного значения функции и. до стороны симплекса определяется как:
Неявная схема (21) консервативна и абсолютно устойчива, имеет четвёртый порядок точности по пространству.
При рассмотрении турбулентных режимов течений, в (1) появляется тензор турбулентных напряжений т.., связанный с видом оператора
У
осреднения. В случае если осреднение незначительное (прямое численное моделирование), т. е. все масштабы течения разрешаются из (1), то тогда т ..^0. Введение осреднения по Рейноль-дсу в (1) приводит к написанию дополнительных сложных моделей для т ., а в (1) рассчитываются только средние величины функций Р и V. Прямое моделирование анизотропной крупномасштабной турбулентной структуры течения напрямую из уравнений Навье-Стокса (1), а осреднение ведётся по малому пространственно-временному масштабу, соответствующему изотропной мелкомасштабной турбулентности (Динамика Больших Вихрей). В таком случае тензор т . имеет простую структуру подсеточной вязкости. Возможна комбинация методов, более детально этот вопрос обсуждается в [11,12]. В данной статье рассматривается модель Динамики Больших Вихрей. Применяется произвольный малый оператор осреднения Ь к (1) по времени и пространству с малым масштабом осреднения, связанным с размерами дискретных элементов в 0х[0.^]. Вводя некоммутативный оператор в Ш4:
V = L(V),
(23)
Тогда в Ш4 имеем разложение вида V=Ь~1(У), коммутирующее над дифференциальным оператором как:
дГ. = Ь-\дУ) = д,(Ь-1 (V)) и т. д. для всех х, (24)
для краткости д1 = д/дt.
В общем случае т . является несимметричным тензором ранга 2, возникающим при применении оператора (23) к (1). В общем случае
тЬ = ¿(V. L~1(V)) - , и точное определение тен-
у J г ^ А зора на малом масштабе осреднения выполняется с
применением модифицированной модели Лэре [6]:
тЬ = VV. - V.V. - 0,5( V. д. V. - V. д. V). (25)
Подставляя (25) в (1) полностью определяем систему уравнений. Легко показать, что применение оператора осреднения по Рейнольдсу
_ л т _
(У = — \VcLt; V' = V - V) вместо (24) в (1) приведёт
к классической незамкнутой цепочке уравнений
7 7
Фридмана-Келлера, в которых т^ = -V) V- .
Для задания корректной начально-краевой задачи для турбулентного течения необходимо учитывать (25). Для большинства задач достаточно во всей области поставить начальный условия (3).
5. Решение уравнения Пуассона на неструктурированной сетке.
Модель учёта свободной поверхности
Для нахождения давления на шаге 3 (2.3) схемы расщепления применяется метод, аналогичный методу вычисления вязких потоков [9]. Окончательно для итерационного нахождения решения уравнения Пуассона приведённая схема меняется на:
1
рп+1 = Щ
i к
akpk
U=1 d/c
-Si
(26)
1
Wt
X
\k=l dk J
где n - номер итерации; K - количество сторон в дискретном симплексе разбиения Q.
Граничные условия типа Неймана реализуются в данной схеме достаточно просто путём подстановки значения производной в (26) саму схему, а граничные условия Дирихле - путём установки соседних "пустых" ячеек.
Наиболее эффективной схемой для решения приведённого уравнения является многосеточный метод. Геометрический многосеточный метод, разработанный автором, описан в [8]. Его плюсы -высокая скорость, быстрая сходимость и малые затраты машинного времени. Минусы - необходимость построения для одной и той же области расчёта набора сеток, в котором идет процедура интерполяции. Алгебраический многосеточный метод находится в стадии разработки.
Для моделирования свободной поверхности используется метод Объёма Жидкостей (Volume of Fluids). Сила поверхностного натяжения определяется как:
F = ok VF, (27)
где o - коэффициент смачивания; k - кривизна свободной поверхности; F - цветовая функция, определяемая уравнением адвекционного переноса в Q:
^ = + к = Щг (28)
Более детально MeTOAVOF рассмотрен в [13].
6. Решение тестовых задач
Необходимо оценить точность метода, область применимости а также численные аспекты предложенной схемы расчёта. Все расчёты выполнены на рабочей станции Лаб.11-3 ИСА РАН: Intel Xeon Quarto 2,33ГГц X2, 8GB RAM.
Для этого рассмотрено решение классической задачи обтекания цилиндра в трёхмерной области. Для расчёта строится адаптивная сетка в области расчёта, достаточная, для разрешения всех важных свойств течения, размер сетки 301 000 элемент. Около поверхности цилиндра используются призмы с целью корректного воспроизведения пограничного слоя. Для расчёта уравнения Пуассона (29) многосеточным методом было построено еще пять более грубых сеток размерами (10, 30, 90, 120, 240>103 элементов. Вначале была проведена оценка скорости сходимости и точности для стационарного режима обтекания цилиндра. Сходимость метода рассмотрена как относительная ошибка по вектор-функции скорости в Lœ, показана на рис.1. Здесь явно видно зависимость от количества элементов используемой сетки.
Для проверки корректности моделирования диффузионной части уравнения рассмотрено стационарное течение около цилиндра для R = 10; 20; 40, Fr = 0. Результаты коэффициентов сопротивления и давления цилиндра сопоставлены с другими публикациями как численных расчётов так и экспериментальных исследований. Сопоставление сведено в таблицу 1. Очевидно, что предложенный численный метод корректно решает данную задачу для малых чисел Рейнольдса.
Для моделирования турбулентного обтекания цилиндра вводится модель Динамики Больших
0,005
* 0,004 ?
N
£ 0,003 ¡ь
I
5 0,002 ■U
Я о.о°1
0,000 О 100 200 300 400
TimeSteps (Г)
Рис.1. Сходимость численного метода в зависимости от количества шагов по времени и количества элементов в сетке
Вихрей (LES), аналогичная (26). Более детально о выборе масштабов осреднения для модели LES можно узнать в [9]. Для турбулентного режима течения выбраны числа Рейнольдса R=(0,5, 1, 10, 15, 25, 55, 85, 100, 110, 350, 1000>104. Выбор обуславливается исследованием кризиса обтекания, который не может быть воспроизведен в случае если численная схема решения НКЗ имеет большую схемную вязкость, поскольку эффекты кризиса обтекания связаны с резкой сметы срыва пограничного слоя в турбулентном режиме около стенки цилиндра (или сферы). Значительная схемная вязкость не позволяет развиться турбулентному пограничному слою. Осреднённые по оси Z значения CD как функция числа Рейнольдса показаны на рис. 2. Экспериментальные данные взяты из [14]. Из рис. 2. явно видно, что предложенный метод обладает возможностью точно разрешать сложные турбулентные течения, чётко описывая режим кризиса обтекания и его последующее восстановление при R > 1-107.
Значения коэффициентов лобового сопротивления Св и коэффициентов давления Ср(х) (я - лоб) для стационарного режима обтекания
Источник Я=10 R=20 Д=40
CD QK0) Ср(п) CD Ср(О) Ср(п) Со ОКО) Ср(п)
[15] 2,828 -0,692 1,500 2,053 -0,582 1,274 1,550 -0,554 1,117
[16] 3,170 -0,687 1,393 2,152 -0,567 1,233 1,499 -0,487 1,133
[17] 3,049 -0,661 1,467 2,048 -0,512 1,289 1,475 -0,448 1,168
Текущий 3,103 -0,666 1,469 2,161 -0,515 1,241 1,505 -0,485 1,143
<Cd>
Рис. 2. Осредненные по высоте (2) значения коэффициента лобового сопротивления Св. Линия - эксперимент [14], восстановленная с бумажной копии работы, X - текущий численный расчёт
Рис. 3. Задача наполнения куба жидкостью с двух сторон изолинии свободной поверхности (цветовой функции)
Для определения точности расчёта задач со свободной поверхностью проведен расчёт модельной задачи наполнения куба с двух сторон жидкостью. В результате расчёта производится проверка закона сохранения массы. Используя гидравлические формулы можно подсчитать объём, попавший в куб к данному моменту времени, и объём, находящийся в кубе в данный момент времени. Расчёт проводился для Ег = 0,1, 0,5, 1,0, 1,5, 3,0, 100,0 при Я = 100; 1-107. Результаты расчёта для Ег = 3,0 и R = 1 • 107 показаны на рис. 3. Для расчёта исходная задача покрывалась 325 000 элементов (все тетраэдры). В результате расчёта на 8 560 шагах по времени (полное наполнение куба) дебаланс массы 0,0015 %, что более чем достаточно при моделировании технических задач.
Кроме показанных задач автором проведены обширные исследования предлагаемого в данной статье численного метода на различных НКЗ.
Построен консервативный и точный метод интегрирования уравнений вязкой несжимаемой жидкости на неструктурированной сетке. Показана высокая точность метода, его безусловная устойчивость и консервативность. Результаты расчёта задач свидетельствуют о проведённых оценках. Алгоритм численного решения систем нелинейный алгебраических уравнений разработан и перенесён на графические вычислительные процессоры. Материалы готовятся к публикации.
Работа поддержана Российским Фондом Фундаментальных Исследований (гранты 08-07-00074а и 09-07-00078а) и программой ОНИТ РАН (проект 1.9).
СПИСОК ЛИТЕРАТУРЫ
1. Белоцерковский О.М. Численное моделирование в механике сплошных сред. М.: Физматлит, 1994.
2. Enright Douglas, Losasso Frank, Fedkiw Ronald.
A Fast and Accurate Semi-Lagrangian Particle Level Set Method // Proceedings of the 4th ASME-JSME Joint Fluids Engineering Conference, № FEDSM2003, 45144. ASME. 2003.
3. Chorda R, Blasco J.A., Fueyo N. An efficient particle-locating algorithm for application in arbitrary 2D and 3D grids // Int. J. of Multiphase Flow. 2002. Vol. 28. № 9. P. 1565-1580.
4. Волков К.Н., Емельянов В.Н. Реализация Ла-гранжевого подхода к описанию течений газа с части-
цами на неструктурированных сетках // Вычислительные методы и программирование. 2008. Т 9. C. 19-33.
5. Paoliy R., Poinsotz T., Shari K. Testing semi-Lagrangian schemes for two-phase flow applications // Proceedings of the Summer Program. Center for Turbulence Research. Toulouse. France. 2006. P. 213-222.
6. Chunlei Liang, Evstigneev N. A study of kinetic energy conserving scheme using finite volume collocated grid for LES of a channel flow // Proceedings of the international conference on numerical methods in fluid dynamics. King's College London, Strand, WC2R 2LS. 2006. P. 61-79.
7. Евстигнеев Н.М., Магницкий Н.А., Сидоров С.В. О природе турбулентности в задаче движения жидкости за уступом // Дифференциальные уравнения. 2009. Т. 45. С. 69-73.
8. Евстигнеев Н.М. Интегрирование уравнения Пуассона с использованием графического процессора технологии CUDA // Вычислительные методы и программирование. 2009. Т. 10. C. 268-274,
9. Евстигнеев Н.М., Магницкий.Н.А., Сидоров С.В. Новый подход к объяснению природы турбулентности вязкой несжимаемой жидкости // Труды ИСА РАН, 2008. Т. 33. C. 49-65.
10. Евстигнеев Н.М. Интегрирование трехмерных уравнений невязкого газа на неструктурированной сетке с применением распределенных вычислений // Вычислительные методы и программирование. 2007. Т. 8. C. 252-264.
11. Cignoni P., Montani C., Scopigno R., Dewall. A
fast divide & conquer Delaunay triangulation algorithm in Ed // Computer J. 2006. Vol. 19. № 2. P. 178-181.
12. Dupont T.F., Liu Y. Back and forth error compensation and correction methods for removing errors induced by uneven gradients of the level set function //
Journal of Computational Physics. 2003. Vol. 190. № 1. P. 311-324.
13. Evstigneev N. Integration of 3D incompressible free surface Navier-Stokes equations on unstructured tetrahedral grid using distributed computation on TCP/IP networks // Proc. Of the VII International conf. "Advances in Fluid Mechanics". Oxford. 2008. 15-20 may. P. 194-208.
14. Roshko A. Experiments on the flow past a circular cylinder at very high Reynolds number // Journal of Fluid Mechanics. 1961. № 10. P. 345-356.
15. Nieuwstadt F. and Keller H.B. Comput. Fluids. 1973. Vol. 1. № 59.
16. He X. and Luo L-S. Phys. Rev. E 55. R6333. 1997.
17. Norberg C. Fluctuating lift on a circular cylinder: Review and new measurements, Journal of Fluids and Structures. 29. 2003. № 17 (1). P. 57-96
18. Evstigneev N.M., Gugushvili I.V. Semi-Lagrangian Method for Advection Equation on GPU in Unstructured Mesh Fluid Dynamics Application // Proc. Of the ICCFD 2009: International Conference on Computational Fluid Dynamics, Bangkok, Thailand during December 25-27. 2009. Vol. 60. P. 1-6.
