Применение метода неполной крестовой аппроксимации к решению задач аэродинамики методом дискретных вихрей Текст научной статьи по специальности «Математика»

УДК 532.5 + 517.9

ПРИМЕНЕНИЕ МЕТОДА НЕПОЛНОЙ КРЕСТОВОЙ АППРОКСИМАЦИИ К РЕШЕНИЮ ЗАДАЧ АЭРОДИНАМИКИ МЕТОДОМ ДИСКРЕТНЫХ ВИХРЕЙ

С.Л. СТАВЦЕВ

Статья представлена доктором технических наук, профессором Ципенко В.Г.1

В статье представлена адаптация метода неполной крестовой аппроксимации для решения нестационарных задач аэродинамики. Новые алгоритмы позволили существенно ускорить расчеты, что подтверждается результатами численных экспериментов на примере задачи трехмерного отрывного обтекания объекта потоком жидкости.

Ключевые слова: несжимаемая жидкость, вихревые методы, нелинейные аппроксимации матриц.

Введение

Вихревые методы являются одной из групп численных методов для решения аэрогидродинамических задач. Одним из наиболее активно развивающихся и широко применяемых направлений среди вихревых методов являются метод дискретных вихрей (МДВ) [1] и его модификации [2]. Метод дискретных вихрей нашел применение в аэродинамике летательных аппаратов, парашютов, в исследованиях ветровой обстановки в условиях городской застройки, при моделировании перемещения смерчей, торнадо.

Ряд решаемых прикладных задач требует значительных вычислительных ресурсов: используемой памяти, времени вычислений. При моделировании вихревых течений методом дискретных вихрей основные вычислительные затраты приходятся на преобразование формы вихревых структур для каждого временного шага. Таким образом, возникает вопрос об ускорении вычислений в вихревых методах.

Ряд методов, направленных на ускорение вычислений скоростей индуцированных вихревыми элементами (вихревой рамкой) в некоторой точке, представлен в работе [3]. В ней повышение эффективности вычислений осуществлялось за счет оправданного снижения точности функции вычисления скорости от каждого вихревого элемента в каждой расчетной точке. Существенным недостатком такого подхода является невозможность контроля точности вычислений, что приводит либо к экспериментальному подбору параметров, влияющих на скорость расчетов, либо к выполнению избыточного числа операций, что сказывается на скорости вычислений.

При решении задач с помощью МДВ более эффективными оказываются методы, использующие иерархические представления данных, например, работы [4; 5]. Основанные на таких идеях методы, которые получили общее название мультипольных методов (fast multipole method - FMM), в последнее время широко применяются для расчета динамики вихревых структур [6; 7]. Преимущество мультипольных методов осуществляется за счет приближенного (часто с заданной точностью) расчета скорости от группы вихревых элементов в нескольких расчетных точках, которые расположены "далеко" друг от друга. Скорости от вихрей, расположенных близко к расчетным точкам, вычисляются явно. Если вихри группировать иерархично,

1 Работа выполнена при поддержке РФФИ, гранты № 11-01-00549-а, 11-01-12137-офи-м-2011, 12-01-00546-а, 10-01-00757-а, гранта Президента МК-6590.2012.1, а также ФЦП "Кадры" Министерства образования и науки, проекты: № 14.740.11.0345, 16.740.12.0727, 14.740.11.1067.

2

то для вычисления скоростей от N вихрей в N расчетных точках вместо 0(К ) потребуется только 0(Мо§^ операций. Несмотря на эффективность мультипольных методов, вычисление мультипольных коэффициентов для каждого класса задач - это нетривиальная подзадача.

Для построения быстрых алгоритмов решения задач аэродинамики эффективным оказывается предложенный в работе [8] иерархический метод мозаично-скелетонной аппроксимации, основанный на вычислении подматриц максимального объема. Этот метод обладает рядом преимуществ перед остальными. Во-первых, для его реализации достаточно знать информацию о взаимном расположении вихрей относительно расчетных точек и функцию вычисления скорости. После этого он может работать как "черный ящик". Во-вторых, для построения приближения достаточно вычислить малое число значений приближаемой функции. В-третьих, приближение строится с заранее заданной точностью.

Алгоритм неполной крестовой аппроксимации был применен в явном виде к решению задач аэродинамики в работе [9]. Эффективность метода можно существенно повысить, если его применить не только для групп вихревых элементов, расположенных близко в пространстве, но и при малых изменениях вихревых элементов и расчетных точек с течением времени, как это сделано в работе [10]. Действительно, в МДВ расчет пелены требуется проводить на каждом временном шаге, поэтому нет необходимости заново строить приближение (аппроксимацию) скорости для удаленных друг от друга групп вихревых элементов и расчетных точек. Более эффективным является использовать данные с предыдущего временного шага. Это может существенно ускорить расчеты.

1. Постановка задачи и алгоритм неполной крестовой аппроксимации

В методе дискретных вихрей [11 - 14] самой трудоемкой операцией является вычисление полей скоростей от вихревых рамок, расположенных на теле, и вихревых элементов пелены в точках вихревой пелены

п

1) = w ¥ (X1) + X — -(х 1),1 = 1,...,ш, (1)

-=1

где w¥(х 1) - невозмущенная скорость набегающего потока; Г - интенсивность --го вихревого элемента (циркуляция рамки или интенсивность вихревого отрезка в пелене); w-(х 1) - скорость, индуцируемая в рассматриваемой точке пелены х 1 вихревым отрезком единичной интенсивности, которая вычисляется согласно закону Био-Савара по формуле

1 г ё1у х(х 1 - У)

— 1) = - 4Рг 1 ГУ-Г". (2)

4Р[х-,х+ ] Iх 1 - У|

В (2) | - криволинейный интеграл по отрезку [х-, х + ] (в (1) четырехугольной рамке [х-,х+ ]

соответствуют 4 слагаемых формулы (2)).

Явный алгоритм расчета поля скоростей (1) выполняется за 0(шп) операций. Для расчета скоростей за 0(шп) операций применим алгоритм неполной крестовой аппроксимации [8; 15]. Сформулируем задачу (1) на матричном языке. На первый взгляд можно найти приближение к матрице, элементы которой являются векторами и вычисляются согласно формуле (2). Однако как показывают численные реализации быстрых алгоритмов, на аппроксимацию матрицы уходит 95 % времени и только 5 % - на умножение полученного приближения на вектор. Поэтому экономнее сжать одну матрицу, даже если это приведет к нескольким умножениям матрицы на вектор.

Также воспользуемся исследованиями работы [3] и на больших расстояниях упростим вычисление скорости (2). В результате построим приближение матрицы А = (а-) с элементами

4p

x i - x j

x i - x j

> a

4plj*(xi -j3

lj •(x i- x +) lj •(x i- x

+

x - x x - x ■

V ij ij

x -x- < a l.

ij j

x + - x ;

^+ - x - x Xj Xj^j

(3)

Параметр а определяется согласно точности аппроксимации. В расчетах параметр выбирался а = 12, что позволяет выполнять вычисления с точностью 10 .

Каждый вектор скорости (1) с помощью матрицы А рассчитывается по формуле

,) = W ¥ (X;) + С; X X; + й;, С, = (Ая); , й ; = (АЬ );, (4)

где а и Ь - вектор-столбцы длиной п с векторными элементами аj j и Ьj = аj X х

] = 1, к, п соответственно.

На первом временном шаге расчет аппроксимации матрицы выполняется согласно стандартному алгоритму, описанному в [8; 15]. При этом вся матрица разбивается на блоки. Если блок соответствует разделенным источникам и приемникам, то его можно сжать, не вычисляя всех его элементов. Процедура сжатия блоков в формализованном виде приведена в виде алгоритма 1.

Алгоритм 1. Построение малоранговой аппроксимации матрицы с заданной точностью. На входе: функция вычисления элементов матрицы А = (а^) (согласно формуле (3)), размер

матрицы т X п, точность вычисляемого приближения е .

A - A

< e AL; U и

На выходе: множители, где A = UV - приближения матрицы A, т.е.

V имеют размеры m х r и n х r соответственно (r - ранг построенного приближения).

p := 1; j1 := |~n/2~|; A := 0 {инициализация: p - номер итерации; j1 - начальный столбец} repeat

u

ip

P-1

:= a-- - Vuikv- k , i = 1,...,m {вычисляем столбец матрицы, по нему очередной вектор

ijp ^^ ik jpk

Up =(Uip ) в U}

k=1

ip = max

U:

vjp := uipj

{определяем номер очередной строки ip}

p-1 >

lipj-1 ui>v

V

p^jk k=1 p у

, j = 1,...,n {вычисляем строку матрицы, по ней очередной век-

= (Vip ) в V}

{определяем номер столбца jp+1 для следующей итерации}

Up =(up-1 up); Vp =(vp-1 Vp) p := p +1

тор vp = vip в Jp+1 = max

il (jmin(m, n) - (p -1) (

until

return U := Up, V := Vp .

U VT

pp

up-1vpt-1 j

> e

U VT

pp

)aND(p < min(m, n))

Из анализа алгоритма 1 видно, что для построения приближения достаточно знать функцию вычисления матричного элемента а-. Таким образом, алгоритм работает как "черный ящик",

что дает ему несомненное преимущество перед остальными методами. Более того, для построения малоранговой аппроксимации нет необходимости вычислять все тп элементов матрицы, а

1

l

3

l

a

1

F

достаточно вычислить только их малую часть, а именно, г(ш + п) элементов. При г << т и г << п значительный выигрыш получается в вычислении вектора у = иУтх : требуется только г(ш + п) действий вместо шп.

2. Модификация алгоритма неполной крестовой аппроксимации

Рассмотрим модификацию алгоритма 1 применительно к задаче динамики. При расчете вихревой пелены расчет скоростей по формуле (1) необходимо выполнять многократно: для каждого временного шага. Однако в прикладных расчетах временные шаги малы и поэтому добавка к аппроксимируемой матрице также будет малой. В этом случае можно не строить малоранговое приближение по алгоритму 1 заново, а использовать аппроксимации с предыдущих

шагов. При этом новое приближение ЛХ+ДХ = ЛХ + ЦУт вычисляется с гораздо меньшими дополнительными рангами, что приводит к значительному уменьшению временных затрат.

Итак, предположим, что на некотором временном шаге X мы построили малоранговое приближение согласно алгоритму 1. На следующем временном шаге X + ДХ необходимо рассчитать скорости

п+Дп

™(х ¡) = ^^ ¥ (х;) + X Г ^ ](Х ¡)Д = 1,...,ш + Дш, (5)

]=1

где Дп - число новых вихревых отрезков; Дш - число новых расчетных точек. В модифицируемом алгоритме выбор новых вихревых отрезков и расчетных точек играет существенную роль: от их выбора зависит ранг дополнения, а значит, скорость работы алгоритма.

Выберем новые расчетные точки и вихревые отрезки для последующего временного шага так, как показано на рис. 1 - в конце вихревой пелены. В этом случае положение остальных отрезков вихревой пелены меняется слабо, поэтому ранг дополнения будет малым.

МОМЫ П" 13 ПМШ Ш /

МОМШ11 ШШМШ1И I ■+ Д/

новые расчетные точки

Рис. 1

В формуле (5) вихревые отрезки и расчетные точки нумеруем таким образом, чтобы номера новых отрезков и точек имели максимальные значения. Тогда матрица ЛХ+ДХ для расчета поля скоростей (5) имеет вид

Л

г+Дг

Л ^

А А

V 21 22 У

где Лг+Д - матрица влияния старых вихревых отрезков на старые расчетные точки. Блоки А12, А21 и Л22 определяются новыми элементами шага и полагаются плотными, т.е. вычисляются явно. Для расчета Лг+Д используем матрицу с предыдущего шага Лг, которая имеет тот же размер. Структуру блоков сжимаемой матрицы Лг+Д1 оставляем без изменений.

Пусть Л0 - некоторый блок матрицы Лг, Лй - соответствующий блок матрицы Л

t+Dt •

Расположение этих блоков и их размеры совпадают. Если блок исходной матрицы A0 плотный, то элементы блока Ad вычисляются явно. Если блок A0 - малоранговый с рангом r0, то из-

T

вестны множители его представления A = U V . В этом случае блок Ad также предполагается малоранговым. Задача заключается в нахождении разложения блока Ad = UV . Так как временной шаг задачи динамики мал, то будем искать решение в виде

U = (и0 Ud); V = (v0 Vd), (6)

где размеры подматриц Ud и Vd m х rd и n х rd соответственно. Для расчета Ud и Vd применяется алгоритм 2, который является модификацией алгоритма 1.

Алгоритм 2. Построение малоранговой аппроксимации дополнения к матрице с заданной точностью.

A0 = U0V0T

На входе: малоранговое представление матрицы A = U V ранга r0; функция вычисления элементов новой матрицы a^; размер матрицы m х n; точность вычисляемого приближения е.

0T

На выходе: матрицы Ud и Vd из (6) такие, что Ad = A + UdVd , rd - ранг приближения

дополнения UdVd .

Р := 1; ji := |~n/2~|; Ad := A0 repeat

Р-1 r0

'Suikvi k -Suikv0k , i = 1,K,m {вычисляем столбец матрицы, по нему очеред

k=1 jp k=1 jp

= (Uip ) в Ud}

uip:= aijp

ной вектор u

ip = max

p i

{определяем номер очередной строки ip }

vjp := uipj

p-1

4pj - Sui"v

k=1

k jk

S u?

0 v0 kvjk

редной вектор vp = vip в Vd

Jp+1

max

i

v

ip

, j = 1,...,n {вычисляем строку матрицы, по ней оче-

k=1 p ' ;

(vip) в Vd}

{определяем номер столбца jp+1 для следующей итерации}

up = (uP-1 up); Vp =(vp-1 vp) p := p +1 until

(Vmin(m,n) - (r0 + p -1) (

U VT

p p

Up-1vpt-1 >e

a0 + UpVpT

)aND(p + r0 < min(m, n))

return Ud := ^ Vd := Vp.

Эффективность алгоритма 2 по сравнению с повторным использованием алгоритма 1 обеспечивается малостью ранга гй по сравнению с рангом г. В результате при построении прибли-

r

жения реже вызывается функция вычисления элемента матрицы, а также выполняется меньшее число вспомогательных операций при построении факторов разложения.

При расчете блоков матрицы на последующих временных шагах в качестве начального приближения используется разложение с множителями (6). Если в результате работы алгоритма на каком-то шаге ранг блока стал + г0 > шт(ш, п), то в дальнейших расчетах этот блок считается плотным. С ростом временных шагов суммарный ранг разреженных блоков накапливается, а также растет число плотных блоков. Это объясняется тем, что из-за смещения вихревых элементов и расчетных точек неизменяемое разбиение матрицы на блоки перестает соответствовать определению малоранговых блоков. Поэтому при выполнении нескольких временных шагов разбиение матрицы на блоки строится заново. Критерий пересчета определяется экспериментальным путем.

3. Построение начального приближения

Одним из ускорений алгоритма 2 является построение хорошего начального приближения к первым г0 -м столбцам искомого малорангового приближения по исходному приближению

А0 = И0 V0 и вычислению малого числа элементов о(г02) аппроксимируемой матрицы. Для нахождения приближения используется метод, изложенный в работе [17], где метод используется для построения приближения трехмерных массивов.

Пусть заданы и0 и V0 , тогда начальное приближение будем искать в виде

А0 = И^0\ (7)

При этом предполагается, что И0 и V0 уже являются хорошим приближением для строк и столбцов блока Ай . Осталось найти матрицу В0. Для ее расчета достаточно вычислить подматрицу размера г0 х г0. В качестве такой подматрицы берется матрица А0, стоящая на пересечении строк и столбцов, отвечающих подматрицам максимального объема И0 и V в исходных матрицах И0 и V0 соответственно. После вычисления А0 из

А ф = И ◊ ВХ (8)

вычисляем В0, обратив И0 и V.

Решение задачи о нахождении подматрицы максимального объема за разумное время представляется затруднительным. В этом случае ищется не наилучшая, а "хорошая" подматрица. В качестве "хорошей" подматрицы И0 матрицы И выбрана такая подматрица, что все элементы

матрицы ИИ-1 не превосходят по модулю 1. Алгоритм нахождения такой подматрицы описан в работе [16].

После вычисления В0 вычисляется ее сингулярное разложение: В0 = ИВX^В (матрицы

ИВ и V с ортонормированными столбцами, XВ - диагональная матрица с ее сингулярными

числами), а также отбрасываются векторы, соответствующие малым сингулярным числам, тогда (7) принимает вид

А0 = И0ИВХВ%^0Т = И~Т , (9)

где И = И0ИВXВ, V = V0VB.

Использование начального приближения в виде (9) позволяет значительно сократить число итераций в алгоритме 2. В алгоритме 2 приближение (9) строится на нулевом шаге.

4. Примеры расчетов

В качестве тестового примера рассматривалось обтекание тонкой пластины под углом атаки 15°. Число ячеек на сетке составляет 9000. Точность расчетов е = 10-3. Результаты расчетов с использованием алгоритмов 1 и 2 представлены на рис. 2, 3. В расчетах первых 10 шагов алгоритмы ускорения расчетов не использовались.

Рис. 2

Рис. 3

На рис. 2 представлены графики зависимости безразмерного времени расчета от временного шага. Время обезразмерено ко времени расчета первого шага алгоритма прямого умножения.

На рис. 3 представлены графики коэффициентов ускорения (отношение времени вычислений прямым методом по сравнению с ускоренным алгоритмом) при использовании алгоритмов 1 и 2 с начальным приближением и без него.

Выводы

В статье приведены алгоритмы расчетов динамики вихревых структур, которые позволяют получить существенное ускорение в расчетах в несколько раз. Это подтверждается результатами численных экспериментов.

ЛИТЕРАТУРА

1. Белоцерковский С.М., Лифанов И.К. Численные методы в сингулярных интегральных уравнениях и их применение в аэродинамике, теории упругости, электродинамике. - М.: Наука, 1985.

2. Дынникова Г.Я. Расчет трехмерных течений несжимаемой жидкости на основе дипольного представления завихренности // Доклады академии наук. - 2011. - Т. 437. - № 1. - С. 35 - 38.

3. Попов В.М., Ненашев А.С. Использование формул приближенных вычислений при решении инженерных задач методом дискретных вихрей // Научный Вестник МГТУ ГА, серия Аэромеханика и прочность. - 2008. - № 125. - С. 102 - 105.

4. Barnes J., Hut P. A hierarchical O( N log N) force-calculation algorithm // Nature, v. 324, № 4, 1986. - P. 446-449.

5. Greengard L., Rokhlin V. A fast algorithm for particle simulations // J. Comput. Physics, v. 73, 1987. - P. 325-348.

6. Дынникова Г.Я. Использование быстрого метода решения "задачи N тел" // Журнал вычислительной математики и математической физики. - 2009. - Т. 49. - № 8. - С. 1458 - 1465.

7. Корнев Н. В. Метод вихревых частиц и его приложение к задачам гидроаэродинамики корабля: дисс. д-ра техн. наук. - СПб.: 1998.

8. Tyrtyshnikov E.E. Mosaic-skeleton approximations // Calcolo, V. 33, 1996. - P. 47-58.

9. Апаринов А.А., Сетуха А.В. О применении метода мозаично-скелетонных аппроксимаций при моделировании трехмерных вихревых течений вихревыми отрезками // Журнал вычислительной математики и математической физики. - 2010. - Т. 50. - № 5. - С. 937 - 948.

10. Stavtsev S.L. Application the method incomplete cross approximation for solving nonstationary dynamics problem of vortexes rings. // Russian Journal of Numerical Analyses and Mathematical Modelling, № 3, 2012 (в печати).

11. Белоцерковский С.М., Ништ М.И. Отрывное и безотрывное обтекание тонких крыльев идеальной жидкостью. - М.: Наука, 1978.

12. Апаринов В.А., Дворак А.В. Метод дискретных вихрей с замкнутыми вихревыми рамками // Труды ВВИА им. проф. Н.Е. Жуковского. - 1986. - Вып. 1313. - С. 424 - 432.

13. Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент. - М.: ТОО "Янус", 1995.

14. Кирякин В.Ю., Сетуха А.В. О сходимости вихревого численного метода решения трехмерных уравнений Эйлера в лагранжевых координатах // Дифференциальные уравнения. - 2007. - Т. 43. - № 9. - С. 1263 - 1276.

15. Tyrtyshnikov E. E. Incomplete cross approximation in the mosaic skeleton method // Computing, v. 64, № 4, 2000. - P. 367-380.

16. Goreinov S.A., Oseledets I.V., Savostyanov D.V., Tyrtyshnikov E.E. and Zamarashkin N.L. How to find a good submatrix. // Matrix Methods: Theory, Algorithms, Applications in V. Olshevsky and E. Tyrtyshnikov editors. - 2010. - P. 247-256.

17. Oseledets I.V., Savostianov D.V., Tyrtyshnikov E.E. Tucker dimensionality reduction of three-dimensional arrays in linear time // SIAM J. Matrix Anal. Appl., V. 30, № 3, 2008. - P. 939-956.

APPLICATION THE METHOD INCOMPLETE CROSS APPROXIMATION FOR SOLVING AERODYNAMIC PROBLEMS BY DISCRETE VORTEXES METHOD

Stavtsev S.L.

In the paper the adaptation of the method incomplete cross approximation for solving nonstationary aerodynamic problems is described. The new algorithms speed up calculation. It is demonstrated by the example of the tree dimensional flow past object.

Key words: vortex methods, nonlinear approximation matrix.

Сведения об авторе

Ставцев Станислав Леонидович, 1978 г.р., окончил ОГУ (2000), кандидат физико-математических наук, старший научный сотрудник ИВМ РАН, автор 27 научных работ, область научных интересов -численные методы решения гиперсингулярных интегральных уравнений, математические методы моделирования задач электродинамики, акустики, аэродинамики.