CC BY
33
Об оценках вычислительной сложности и погрешности быстрого алгоритма в методе вихревых элементов
К.С. Кузьмина <[email protected]> И.К. Марчевский <[email protected] > МГТУ им. Н.Э. Баумана, 105005, Россия, г. Москва,ул. 2-я Бауманская, д. 5
Аннотация. Основная вычислительная сложность при использовании вихревых методов сосредоточена в вычислении конвективных и диффузионных скоростей всех вихревых элементов. Рассматривается один из эффективных способов ускорения вычислений в методе вихревых элементов — алгоритм типа Барнса - Хата. Метод основан на построении иерархической структуры областей (дерева). При вычислении конвективных скоростей вихревых элементов данный метод позволяет учитывать взаимное влияние вихревых элементов, находящихся далеко друг от друга, приближенно по формуле, полученной с помощью разложения выражения для конвективной скорости в ряд Тейлора. Влияние вихревых элементов, находящихся в ближней зоне, рассчитывается напрямую по закону Био — Савара. При реализации данного алгоритма возникают параметры, влияющие на вычислительную трудоемкость и точность решения задачи: к — количество уровней дерева и в — параметр дальности ячеек. Влияние вихревых элементов на диффузионные скорости друг друга экспоненциально затухает с увеличением расстояния между ними. Поэтому для вычисления диффузионных скоростей также построен алгоритм, позволяющий с помощью использования структуры дерева находить вихревые элементы, находящиеся в ближней зоне и вычислять влияние только от них. На основе решения модельных задач получены оценки вычислительной сложности алгоритмов вычисления конвективных и диффузионных скоростей, которые зависят от параметров алгоритма и количества вихревых элементов. Также получены оценки погрешности вычисления конвективных и диффузионных скоростей в зависимости от параметров алгоритма. На практике эти оценки позволяют выбирать оптимальные значения параметров алгоритма и добиваться максимального ускорения вычислений при заданном уровне допустимой погрешности расчета.
Ключевые слова: метод вихревых элементов, закон Био - Савара, вязкая жидкость, диффузионная скорость, алгоритм Барнса - Хата, быстрый мультипольный метод, вычислительная сложность, оценка погрешности.
БО!: 10.15514/18РЯЛ8-2016-28(1 )-15
Для цитирования: Кузьмина К.С., Марчевский И.К. Об оценках вычислительной сложности и погрешности быстрого алгоритма в методе вихревых элементов. Труды ИСП РАН, том 28, вып. 1, 2016 г., с. 259-274. DOI: 10.15514/ISPRAS-2016-28(1)-15
1. Введение
Метод вихревых элементов — бессеточный лагранжев метод вычислительной гидродинамики, позволяющий моделировать течения идеальной или вязкой несжимаемой среды. Этот метод весьма эффективен при моделировании внешнего обтекания профилей, когда область течения является неограниченной. Основными достоинствами вихревых методов при решении таких задач являются низкая численная диссипация и точное выполнение граничного условия на бесконечности, в то время как при использовании сеточных методов для моделирования внешнего обтекания необходимо ограничивать расчетную область и обеспечивать удовлетворение некоторых «искусственных» граничных условий. Кроме того, вихревые методы позволяют моделировать не только течения вокруг неподвижных жестких тел, но также и гидроупругие колебания тел, сохраняя при этом примерно ту же вычислительную сложность.
Вихревые методы достаточно популярны в различных инженерных приложениях. Они позволяют вычислять как параметры течения, так и аэрогидродинамические силы, действующие на тело, с достаточной точностью, при этом их вычислительная сложность может быть значительно ниже по сравнению с сеточными методами, особенно при использовании в вихревых методах приближенных быстрых алгоритмов. Целью данной работы является построение оценки вычислительной сложности быстрой реализации алгоритма вихревых методов и исследование зависимости погрешности приближенного решения от времени выполнения расчета.
2. Краткое описание вихревых методов
Существует несколько модификаций вихревых методов для моделирования двумерных и пространственных течений, иногда существенно отличающихся друг от друга как в части реализации отдельных операций и алгоритмов, так и в смысле используемых методов и подходов. Однако их общей особенностью является то, что первичной расчетной величиной является завихренность, а ее распределение в области течения моделируется большим количеством отдельных вихревых элементов. Каждый вихревой элемент генерирует «элементарное» поле скоростей во всей области течения, а общее поле скоростей является суперпозицией этих «элементарных» полей и выражается с помощью закона Био - Савара:
1 ГП(^)х(г-0
':] -^^-^ (2Л)
5 |Г — ^
2(а-1)п.
Здесь 5 — область с ненулевой завихренностью; й — размерность пространства в задаче; £:) = V х £:) — вектор завихренности в точке с радиус-вектором
Способ аппроксимации распределения завихренности определяется выбором модели вихревого элемента. В двумерных задачах обычно используется модель круглого вихря Рэнкина; для пространственных задач известно несколько моделей вихревых элементов, каждая из которых имеет свои достоинства и недостатки. Далее будем рассматривать двумерный случай, но основные идеи после некоторой адаптации могут быть применены и к решению пространственных задач с помощью вихревых методов. Течение вязкой несжимаемой жидкости описывается с помощью уравнений Навье - Стокса, которые могут быть записаны в форме Гельмгольца:
дИ ^
— + (7 • V)n + уДП = 0,
где V — коэффициент кинематической вязкости.
Наиболее простым способом учета конвективного слагаемого
(7 • V)n
является перемещение всех вихревых элементов вдоль линий тока, при этом их интенсивности (циркуляции) остаются постоянными; такой алгоритм будет соответствовать моделированию течения идеальной жидкости. Для учета слагаемого, отвечающего за влияние вязкости, известно несколько подходов: метод случайных блужданий [1], метод обмена интенсивностями между вихревыми элементами [2], метод диффузионной скорости [3,4,5]. В соответствии с подходом, предполагающим введение диффузионной скорости, уравнение Навье - Стокса в двумерном случае может быть записано в следующей форме:
Здесь ^^ = ^ V? — диффузионная скорость. Таким образом, при моделировании вязкого течения методом вязких вихревых доменов (ВВД) [4] необходимо численное решение следующей системы на каждом временном шаге:
'¿Г;
-у- = 0,
^ ¿ = 1.....N.
= + ТО),
Первое уравнение означает, что интенсивности вихревых элементов Г остаются постоянными (как в случае идеальной жидкости), а их положения меняются таким образом, что они движутся вдоль линий тока поля скоростей V + ЙЛ
3. Вычислительная сложность метода вихревых элементов
Основная вычислительная сложность вихревых методов сосредоточена в вычислении конвективных и диффузионных скоростей всех вихревых элементов. Для вихревых элементов Рэнкина согласно закону Био -Савара (2.1) получаем следующее выражение для конвективной скорости
-> ^ -> V Г кх (г — ш
¿=1
Здесь е — радиус вихревого элемента, который полагается постоянным; т^ — положение -го вихревого элемента; Г — его интенсивность; N — количество вихревых элементов. Ясно, что вычислительная сложность процедуры вычисления конвективных скоростей пропорциональна Ы2. Если учитывать только операции умножения и деления, она равна 6Ы2. Можно добиться уменьшения сложности, если учесть, что Ку и имеют противоположные знаки и отличаются только множителями Г, однако сложность в любом случае будет превышать 3Ы2 — вклад -го вихревого элемента в конвективную скорость ¿-го вихревого элемента).
Прямое вычисление скоростей (с использованием формулы (3.1)) возможно, когда количество вихревых элементов не превышает нескольких десятков тысяч, в противном случае время выполнения расчета становится недопустимо большим. Обычно для моделирования нестационарного течения, особенно в гидроупругих задачах, необходимо выполнять до нескольких тысяч шагов по времени; для повышения точности моделирования течения необходимо увеличивать количество вихревых элементов, что в свою очередь приведет к увеличению времени расчета одного шага по времени. В данной работе не затрагивается проблема моделирования обтекания профиля, т.е. не рассматривается учет граничных условий прилипания (непротекания в случае идеальной жидкости) на поверхности профиля. Отметим лишь, что для учета влияния профиля на течение в вихревых методах профиль заменяется вихревым слоем (и слоем источников в случае движущегося или деформируемого профиля), интенсивность которого может быть определена из решения некоторого интегрального уравнения [6,7,8]. В общем случае вихревой слой является свободным, его также можно моделировать с помощью вихревых элементов, которые становятся частью вихревого следа.
Для снижения временных затрат на выполнение вычислений используются разные методы ускорения расчетов. Один из наиболее популярных подходов — использование алгоритма типа Барнса - Хата [9], изначально разработанного для приближенного быстрого решения гравитационной задачи N тел. К вихревым методам этот алгоритм был адаптирован в работе [10]. Этот алгоритм основан на построении иерархической структуры областей — дерева. Напрямую по закону Био - Савара (3.1) вихревое влияние рассчитывается только от вихревых элементов, которые находятся в достаточно близких к рассматриваемой ячейках, влияние от вихревых элементов, находящихся далеко от рассматриваемой ячейки, рассчитывается приближенно по упрощенным формулам.
Используя данный подход, можно добиться значительного снижения вычислительной сложности алгоритма: при правильном выборе количества уровней дерева сложность будет пропорциональна NlogN. Это позволяет за разумное время решать задачи с большим количеством (десятками или даже сотнями тысяч) вихревых элементов.
При решении практических задач недостаточно знать только порядок вычислительной сложности, необходима более подробная оценка, которая может помочь при выборе оптимальных параметров алгоритма. Известные оценки [11, 12] позволяют выбирать оптимальное количество уровней дерева, но они имеют невысокую точность. В [13] авторами получена более точная оценка количества операций:
* 2к Ы( Vw ) ( e(V2)4 Vw 1 1
896-2fc-£/ / 1 1 \ /(V2)fc-4\\
Здесь N — количество вихревых элементов в области течения, к — количество уровней дерева, в — параметр точности. Значения коэффициентов а и р определяются эмпирически, они зависят от решаемой задачи. Значение а
(0 < а < 1) значительно зависит от распределения вихревых элементов в вихревом следе, а величина р (^ > 0) зависит от формы вихревого следа. При малых значениях параметра в (0 < в << 1) быстрый метод имеет высокую точность, но в этом случае он имеет также и высокую вычислительную сложность; большие значения в позволяют понизить сложность вычислений, но при этом повышается погрешность решения.
Известны работы, в которых получены и математически доказаны оценки погрешности алгоритма Барнса - Хата [14,15], но эти оценки имеют асимптотический характер и содержат некоторые неизвестные параметры. Эти
оценки могут быть полезны для «теоретической» оценки, но на практике их применение затруднительно. Целью данной работы является построение приближенных оценок для вычислительной сложности быстрого алгоритма, а также для погрешности вычисления скоростей.
4. Модельная задача
Рассмотрим тестовую задачу о моделировании течения вязкой жидкости, которая соответствует известному явлению диффузии круглого вихря в неограниченной области (вихря Ламба). Пусть в начальный момент времени £ = 0 в вязкую несжимаемую жидкость введена бесконечная вихревая нить с циркуляцией Г. Точное решение данной задачи для распределения завихренности в момент времени £ имеет вид
где г — расстояние до центра вихря. Найдем суммарную завихренность внутри окружности радиуса Д:
2л: й ,
rfi(t) = J d<p J fi(r, t) r dr = Г (1 — exp ( -
д2у
Будем считать, что Г = 1, коэффициент кинематической вязкости среды V = (2000тс)-1 « 0,00016, время начала моделирования распределения завихренности £0 = 2000^ « 6283. Тогда более 99,8 % всей завихренности содержится в окружности радиуса Д = 5, остальной пренебрежем. Количество вихревых элементов, которое было использовано для моделирования распределения завихренности внутри этого круга, взято из интервала N = 60 000...300 000, при этом их пространственное распределение близко к равномерному. Пример такого распределения (с небольшим количеством вихревых элементов N « 1000) показано на рис. 1. Интенсивности вихревых элементов вычисляются аналитически при помощи интегрирования точного решения (4.1).
Рис.1. Распределение вихревых элементов в круглом вихре (N ~ 1000) Fig. 1. The distribution of the vortex elements in a circular vortex
0
0
Для построения дерева для алгоритма типа Барнса - Хата был использован стандартный метод [10], ячейки дерева для разных уровней дерева показаны на рис. 2.
2-й уровень 4-й уровень
6-й уровень 8-й уровень
Рис. 2. Ячейки для разных уровней дерева Fig.2. Cells for different levels of the tree
5. Оценка вычислительной сложности
Для построения эффективного алгоритма необходима подробная оценка вычислительной сложности всех его частей. Отметим, что построение дерева является очень «легкой» процедурой (даже при количестве вихревых элементов порядка сотен тысяч) по сравнению с вычислением конвективных и диффузионных скоростей всех вихревых элементов.
5.1. Вычисление конвективных скоростей
Оценка вычислительной сложности алгоритма расчета конвективных скоростей состоит из двух частей: первая соответствует прямому вычислению скоростей с помощью закона Био - Савара от вихревых элементов, ячейки которых расположены достаточно близко друг к другу, согласно [13] эта оценка имеет следующую форму:
24N2 /4\:2 (V2) - 1
2
1 —-
0(V2)K
'l — ^1
Vw
+ 4N.
(5.1)
Вторая часть оценки касается приближенного вычисления скоростей от ячеек, которые находятся достаточно далеко от рассматриваемой ячейки:
896-2fc-S/ / 1 1
^Far =-- ( 4 ( + ■
+ ln|(V^l). (5.2)
02 V 14 + б 4-(^б/ V 4 + б „
Критерий дальности ячеек имеет следующую форму:
в • 5 > й + й0,
где 5 — октаэдрическая норма (У • ||1) вектора, соединяющего центры «влияющей» и «рассматриваемой» ячеек, й и й0 — суммы длин сторон этих ячеек.
Для того чтобы получить оценку коэффициента а в формулах (3.2) и (5.1), с помощью специальной процедуры в алгоритме был произведен точный подсчет количества операций, осуществляемых в расчетах. Оценка для (5.1) содержит только параметр а, а параметр р не входит в это выражение, что позволяет получить хорошую оценку для параметра а. Результаты подсчета количества операций с помощью численного
эксперимента (реальное количество операций) и с помощью аналитической оценки (5.1) для 150 000 вихревых элементов при значении параметра а = 0,8 показаны на рис. 3.
Рис. 3. Количество операций (N = 150 000, а = 0,8) Fig. 3. The number of operations QRS (N = 150000, a = 0,8)
4
Для получения оптимального значения параметра а были проведены численные эксперименты с разным количеством вихревых элементов (М = 60 000 ... 300 000), моделирующих исходный вихрь, для разных значений параметра в = (0,1... 1,0). Количество уровней дерева при этом выбиралось исходя из оценок, приведенных в работах [12, 13]. Для каждого значения N было вычислено значение коэффициента а, минимизирующее сумму квадратов относительных погрешностей:
где б; = 0,1i; — количество операций, полученное из численного
эксперимента для N вихревых элементов и в =0;.
Вычисленные оптимальные значения параметра а приведены в табл. 1. Полученные значения близки друг к другу, поэтому на практике можно использовать среднее значение а = 0,844.
Похожая процедура подсчета количества операций и минимизации погрешностей была проведена для вычисления оптимального значения параметра р в оценке (5.2). Полученные оптимальные значения параметра р для разных N приведены в табл. 1. Снова можно заметить, что эти величины близки, поэтому предлагается, так же как и в случае с а, использовать среднее значение р = 0,561.
Табл. 1. Оптимальные значения а и р для разных значений N (fc — количество уровней дерева)
Table 1. Optimal values a and р for different N (к - the number of levels of the tree)
N 60 000 90 000 120 000 150 000 180 000 210 000 240 000 270 000 300 000
~k 14 15 15 16 16 16 16 16 17
0,885 0,793 0,826 0,837 0,863 0,876 0,871 0,876 0,776
P 0,521 0,558 0,587 0,539 0,554 0,566 0,582 0,587 0,582
Таким образом, получена оценка вычислительной сложности алгоритма вычисления конвективных скоростей.
5.2. Вычисление диффузионных скоростей
Вычисление диффузионных скоростей производится с помощью похожего подхода: во внимание принимается влияние только от тех вихревых элементов, которые находятся в ячейках нижнего уровня, расположенных достаточно близко к рассматриваемой ячейке и удовлетворяющих условию:
0di/ • (5 - 0.5(h + he)) < £*,
где — параметр дальности ячеек для диффузионных скоростей; чем меньше этот параметр, тем точнее приближенное решение; 5, h и h0 — те же величины, что и в (5.3); е* — характерное расстояние между вихревыми элементами. Расчетная формула для вычисления диффузионной скорости приведена, к примеру, в [4].
Для того, чтобы получить выражение для оценки вычислительной сложности алгоритма расчета диффузионных скоростей, можно использовать те же рассуждения, что и для конвективных скоростей [13]. Опуская промежуточные вычисления, отметим, что вычислительная сложность для алгоритма расчета диффузионных скоростей пропорциональна сложности Q алгоритма для конвективных скоростей (3.2), но при этом в нужно заменить на в
di/ •
Y
Qdi/ = @le=edif •N •
Здесь N и fc, как и ранее, количество вихревых элементов и уровней дерева; у — коэффициент, который может быть оценен численно (табл. 2).
Табл. 2. Оптимальные значения коэффициента у для разных значений N (к — количество уровней дерева)
Table 2. Optimal values of coefficient у for different N (k - the number of levels of the tree)
N 60 000 90 000 120 000 150 000 180 000 210 000 240 000 270 000 300 000
k 14 15 15 16 16 16 16 16 17
У 0,0775 0,0775 0,0731 0,0712 0,0694 0,0689 0,0686 0,0679 0,0686
Видно, что значения коэффициента у при различных N близки, поэтому в расчетах будем использовать среднее значение у ~ 0,7.
6. Погрешности вычисления скоростей с помощью быстрых алгоритмов
Основной целью применения алгоритма Барнса - Хата является ускорение вычислений. Полученные формулы (5.1), (5.2) и (5.4) позволяют оценить, сколько времени требуется для моделирования течения, но кроме этого немаловажным вопросом является оценка величины погрешности быстрого метода. Далее будем рассматривать относительную погрешность вычисления конвективных и диффузионных скоростей
£ =
max — 7*
¿=1.....w1 ' '
где — конвективная либо диффузионная скорость -го вихревого элемента, вычисленная быстрым методом; — скорость того же вихревого элемента, вычисленная напрямую по формуле (3.1) для конвективных скоростей и по формуле [4] для диффузионных скоростей; | Ус*опр | — максимальное значение конвективных скоростей всех вихревых элементов. Для модельной задачи, которая рассматривается в данной работе, ~ 0,05 и это значение
практически не зависит от количества вихревых элементов N.
6.1. Погрешность вычисления конвективных скоростей
Погрешность вычисления конвективной скорости значительно зависит от значения параметра в; большие значения параметра в соответствуют большим погрешностям.. Для того, чтобы оценить, как параметр в влияет на погрешность вычисления конвективной скорости, был выполнен ряд вычислительных экспериментов. Были произведены расчеты для N = 60 000...300 000 с использованием быстрого метода с разными параметрами в = 0.1... 1.0. Полученные результаты изображены на рис. 4; тонкие линии соответствуют относительным погрешностям вычисления конвективной скорости для экспериментов с разными М, жирная линия соответствует кривой, которая является их оценкой сверху.
Рис. 4. Относительная погрешность для конвективных скоростей (тонкие линии соответствуют N = 60 000 ... 300 000, жирная линия — оценка сверху) Fig 4. The relative error for the convective velocity (thin lines correspond to N = 60 000 ... 300 000, thick line - the upper bound)
Можно заметить, что для каждого значения N зависимость погрешности от в хорошо аппроксимируется функцией econv = с •б3, при этом для разных N значение коэффициента с варьируется от 0,020 до 0,025. Таким образом, мажоранта этих кривых
econv = 0,025 •б3
соответствует достаточно точной оценке погрешности вычисления конвективных скоростей быстрым методом.
6.2. Погрешность вычисления диффузионных скоростей
Погрешность вычисления диффузионной скорости с помощью быстрого метода зависит не только от параметра , но и от некоторых других параметров. Для каждого N относительная погрешность близка к квадратичной функции со своим коэффициентом пропорциональности. На рис. 5 показаны соответствующие кривые для N = 90 000 и N = 300 000. Соответствующий коэффициент пропорциональности может быть найден с помощью метода наименьших квадратов.
Рис. 5. Относительные погрешности для диффузионных скоростей (тонкие линии —
численный эксперимент; жирные линии — квадратичные функции £ = )
Fig. 5. The relative error for the diffusion rates (thin line - numerical experiment; thick lines - quadratic functions £ =
Для больших значений параметра наблюдается значительное различие между квадратичной оценкой и результатами численных экспериментов, но с практической точки зрения это не имеет большого значения, поскольку погрешность в 10 - 20 % недопустима и необходимо производить расчеты при малых значениях параметра , для которых оценка и результаты эксперимента хорошо согласуются.
Значения коэффициента , вычисленные для разных экспериментов (W = 60 000 ... 300 000), представлены в табл. 3.
Табл. 3. Значения коэффициента для разных значений N (fc — количество уровней дерева)
Table 3. Values of the coefficient for different N (k - the number of levels of the tree)
N 60 000 90 000 120 000 150 000 180 000 210 000 240 000 270 000 300 000
k 14 15 15 16 16 16 16 16 17
cdi/ 6,86 9,81 9,76 13,70 13,76 13,19 13,51 13,72 19,28
Из табл. 3 видно, что значение коэффициента зависит только от
количества уровней дерева, т.е. = а для фиксированного значения
к оно практически постоянно. Можно попробовать записать приближенную формулу для (используя метод наименьших квадратов или какой-либо
другой способ), но она будет чисто эмпирической и сложно применимой на практике. Получение аналитической оценки этого явления является целью дальнейших исследований.
7. Заключение
Описана подробная оценка вычислительной сложности алгоритма типа Барнса - Хата для процедуры вычисления конвективных и диффузионных скоростей в методе вихревых элементов. Получена зависимость погрешности вычисления конвективных и диффузионных скоростей для модельной задачи. Если принять допустимой погрешность е = 0,2 %, то можно выбрать следующие значения параметров дальности для алгоритма быстрого метода: в = 0,4, = 0,1. Посчитаем суммарную вычислительную сложность расчета конвективных и диффузионных скоростей быстрым методом Qtot и сравним ее со сложностью прямого метода = 15N2 (коэффициент 15 получен с учетом того, что расчет одного парного влияния на конвективные скорости требует 6 операций умножения и деления, а на диффузионные — 9 операций). Величина ускорения 5 = показывает эффективность
использования быстрого метода (табл. 4).
Табл. 4. Вычислительные сложности прямого и быстрого метода и ускорение Table. 4. Computational complexity of direct and rapid methods and the acceleration
N 60 000 90 000 120 000 150 000 180 000 210 000 240 000 270 000 300 000
k 14 15 15 16 16 16 16 16 17
109 0,58 0,94 1,26 1,70 1,93 2,24 2,60 2,97 3,69
109 54,3 122,0 215,9 343,2 486,4 670,5 864,3 1097,4 1357,0
s 93,6 129,8 171,3 201,9 252,0 299,3 332,4 369,5 367,8
Таким образом, при N « 70 000 достигается 100 кратное ускорение вычислений по сравнению с прямым методом. А при расчете N « 250 000 можно добиться 350 кратного ускорения.
Полученные оценки позволяют выбирать параметры для алгоритма быстрого метода таким образом, чтобы получить минимально возможную вычислительную сложность для требуемой величины погрешности.
Список литературы
[1]. A.J. Chorin, Numerical study of slightly viscous flow, J. Fluid. Mech., 57 (1973), pp. 785-796.
[2]. P. Degond, and S. Mas-Gallic, The weighted particle method for convection-diffusion equations. Part I: The case of an isotropic viscosity, Math. Comp., 53 (1989), pp. 485507.
[3]. Ogami Y., Akamatsu T. Viscous flow simulation using the discrete vortex model — the diffusion velocity method // Computers and Fluids. 1991. V. 19, N 3/4. P. 433-441.
[4]. Дынникова Г.Я. Лагранжев подход к решению нестационарных уравнений Навье -Стокса // Доклады Академии наук. 2004. Т. 399, №1. С. 42-46.
[5]. S. Guvernyuk, and G. Dynnikova, Modeling the ow past an oscillating airfoil by the method of viscous vortex domains, Fluid Dynamics, 42 (2007), pp. 1-11.
[6]. I. K. Lifanov, and S. M. Belotserkovskii, Methods of Discrete Vortices. CRC Press, 1993.
[7]. S. N. Kempka, M. W. Glass, J. S. Peery, and J. H. Strickland, Accuracy Considerations for Implementing Velocity Boundary Conditions in Vorticity Formulations. SANDIA Report SAND96-0583, 1996.
[8]. K. S. Kuzmina, and I. K. Marchevsky The Modified Numerical Scheme for 2D Flow-Structure Interaction Simulation Using Meshless Vortex Element Method, in Proc. IV Int. Conference on Particle-based Methods — Fundamentals and Applications (PARTICLES-2015), Barcelona (2015), pp. 680-691.
[9]. J. Barnes, and P. Hut, A hierarchical O(WlogW) force-calculation algorithm, Nature, 324 (1986), pp. 446-449.
[10]. Дынникова Г.Я. Использование быстрого метода решения «задачи N тел» при вихревом моделировании течений // Журнал вычислительной математики и математической физики. 2009. Т. 49, № 8. С. 1458-1465.
[11]. Гирча А.И. Быстрый алгоритм решения «задачи N тел» в контексте численного методя вязких вихревых доменов // Информационные технологии моделирования и управления. 2008. № 1. С. 47-52.
[12]. Морева В.С. Способы ускорения вычислений в методе вихревых элементов // Вестн. Моск. гос. техн. ун-та им. Н.Э. Баумана. Сер.: Естественные науки. Спец. выпуск «Прикладная математика». 2011. С. 83-95.
[13]. Кузьмина К.С., Марчевский И.К. Оценка трудоемкости быстрого метода расчета вихревого влияния в методе вихревых элементов // Наука и образование: электрон. науч.-техн. изд. 2013. № 10. С. 399-414.
[14]. A. Grama, V. Sarin, and A. Sameh, Improving Error Bounds for Multipole-Based Treecodes, SIAM J. Sci. Comp., 21 (2000), pp. 1790-1803.
[15]. J. K. Salmon, and M. S. Warren, Skeletons from the treecode closet, J. Comput. Phys., 111 (1994), pp. 136-155.
On the estimations of efficiency and error of fast algorithm in vortex element method
K.S. Kuzmina <[email protected] > I.K. Marchevsky <[email protected] > Bauman Moscow State Technical University, 105005, Russia, Moscow, 2-nd Baumanskaya st., 5.
Abstract. The main computational complexity of vortex methods is concentrated in the calculation of the convective and diffusive velocities of all vortex elements. The analogue of the Barnes - Hut algorithm is considered as one of the most efficient ways to acceleration of the velocities computation in vortex element method. This method is based on the tree (hierarchical structure of regions) construction. When calculating the convective velocities this algorithm makes it possible to take into account the influence of vortex elements, which are located far enough from each other, approximately by using Taylor approximation of expression for convective velocity. The influence of vortex elements, which are close to each other is calculated directly using Biot —Savart law. There are two parameters of algorithm that affect the computational complexity of the algorithm and the problem solving accuracy: fc is the maximum tree level, 0 is the parameter which determines the radius of a near-field zone. When calculating diffusive velocities the influence of the vortex elements to each other decays exponentially with increasing distance between them. Therefore, constructed algorithm of diffusive velocities calculation allows finding vortex elements from near-field zone using tree structure and calculating diffusive velocities using only vortex elements from this zone.
The estimations of computational complexity, which depends on algorithm parameters and number of vortex elements, are obtained for the algorithms for convective and diffusive velocities calculation. Also estimations for the errors of the vortex elements velocities computation, which depend on the algorithm parameters, are constructed. These estimates allow in practice to choice the optimal parameters of the whole algorithm and achieve the maximum possible acceleration of the computations for the given maximum error level.
Keywords. Vortex element method, Biot - Savart law, viscous flow, diffusive velocity, Barnes - Hut algorithm, fast multipole method, computational complexity, error estimation
DOI: 10.15514/ISPRAS-2016-28(1 )-15
For citation: Kuzmina K.S., Marchevsky I.K. On the estimations of efficiency and error of fast algorithm in vortex element method. Trudy ISP RAN/Proc. ISP RAS, 2014, vol. 28, issue 1, pp. 259-274 (in Russian). DOI: 10.15514/ISPRAS-2016-28(1)-15
References
[1]. A.J. Chorin, Numerical study of slightly viscous flow, J. Fluid. Mech., 57 (1973), pp. 785-796.
[2]. P. Degond, and S. Mas-Gallic, The weighted particle method for convection-diffusion equations. Part I: The case of an isotropic viscosity, Math. Comp., 53 (1989), pp. 485507.
[3]. Y. Ogami, and T. Akamatsu, Viscous flow simulation using the discrete vortex model — the diffusion velocity method, Computers & Fluids, 19 (1991), pp. 433-441.
[4]. G. Ya. Dynnikova, Lagrange method for Navier - Stokes equations solving, Doklady Akademii Nauk, 399 (2004), pp. 42-46.
[5]. S. Guvernyuk, and G. Dynnikova, Modeling the flow past an oscillating airfoil by the method of viscous vortex domains, Fluid Dynamics, 42 (2007), pp. 1-11.
[6]. I. K. Lifanov, and S. M. Belotserkovskii, Methods of Discrete Vortices. CRC Press, 1993.
[7]. S. N. Kempka, M. W. Glass, J. S. Peery, and J. H. Strickland, Accuracy Considerations for Implementing Velocity Boundary Conditions in Vorticity Formulations. SANDIA Report SAND96-0583, 1996.
[8]. K. S. Kuzmina, and I. K. Marchevsky, The Modified Numerical Scheme for 2D Flow-Structure Interaction Simulation Using Meshless Vortex Element Method, in Proc. IV Int. Conference on Particle-based Methods — Fundamentals and Applications (PARTICLES-
2015), Barcelona (2015), pp. 680-691.
[9]. J. Barnes, and P. Hut, A hierarchical O(WlogW) force-calculation algorithm, Nature, 324 (1986), pp. 446-449.
[10]. G. Ya. Dynnikova, Fast technique for solving the W-body problem in flow simulation by vortex methods, Computational Mathematics and Mathematical Physics, 49 (2009), pp. 1389-1396.
[11]. A. I. Gircha, Fast Algorithm for W-body Problem Solving with Regard to Numerical Method of Viscous Vortex Domains, Informatial technologies in Simulating and Control, 1 (2008), pp. 47-52.
[12]. V. S. Moreva, On the Ways of Computations Acceleration when Solving 2D Aerodynamic Problems by Using Vortex Element Method, Heralds of the Bauman Moscow State University. Natural Sciences. Sp. Issue 'Applied Mathematics' (2011), pp. 83-95.
[13]. K. S. Kuzmina, I. K. Marchevsky, Estimation of computational complexity of the fast numerical algorithm for calculating vortex influence in the vortex element method, Science & Education (electronic journal), 10 (2013), pp. 399-414 (URL: http://technomag.bmstu.ru/en/doc/604030.html)
[14]. A. Grama, V. Sarin, and A. Sameh, Improving Error Bounds for Multipole-Based Treecodes, SIAM J. Sci. Comp., 21 (2000), pp. 1790-1803.
[15]. J. K. Salmon, and M. S. Warren, Skeletons from the treecode closet, J. Comput. Phys., 111 (1994), pp. 136-155.
