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

Том ХЫ

УЧЕНЫЕ ЗАПИСКИ ЦАГИ 2010

№ 1

УДК 534.2

519.63:517.96

МЕТОДИКА ПОВЫШЕНИЯ ТОЧНОСТИ ПРИ МОДЕЛИРОВАНИИ ПЕРЕНОСА АКУСТИЧЕСКИХ ВОЗМУЩЕНИЙ НА НЕСТРУКТУРИРОВАННЫХ СЕТКАХ

И. В. АБАЛАКИН, А. ДЕРВЬЁ, Т. К. КОЗУБСКАЯ, Х. УВРАР

Рассматриваются два пути повышения точности численного моделирования переноса акустических возмущений на неструктурированных сетках при использовании алгоритмов с определением переменных в узлах сетки. Апостериорная оценка точности обеих методик проводится на примере численного решения модельной задачи.

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

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

С математической точки зрения точность алгоритма при условии гладкости решения полностью определяется точностью устойчивой аппроксимации искомой функции. Условия реальных аэроакустических приложений, характеризуемых трехмерной геометрией, сильным перепадом масштабов, присутствием криволинейных поверхностей и т. п., накладывают свои требования на сетку, что приводит в результате к ее достаточно сложной структуре. Так, в данной работе речь идет о неструктурированных треугольных и тетраэдральных сетках. Обеспечение высокой точности аппроксимации на таких сетках сопряжено с серьезными трудностями и представляет отдельную проблему. Именно этой проблеме посвящена настоящая статья. В ней рассматриваются два подхода, построенные на основе метода конечных объемов, к повышению точности алгоритмов для моделирования распространения акустических возмущений на неструктурированных треугольных и тетраэдральных сетках при определении переменных в узлах сетки*. Первый — это разрабатываемые авторами многопараметрические схемы повышенной точности [1—4]. Базовый вариант схемы использует линейную реконструкцию переменных на границах соседних ячеек [3]. Второй подход основан на внедрении полиномиальной реконструкции [5].

1. Постановка задачи. Так как данная работа посвящена алгоритмам для моделирования переноса акустических возмущений, то в качестве основной математической модели, не ограничивая общности, рассмотрим систему линеаризованных уравнений Эйлера в двумерной декартовой системе координат. В случае стационарного и однородного по пространству среднего течения она имеет вид:

* Условие определения переменных в узлах является принципиальным для дальнейших рассуждений.

дР + (^, У)р' + рdivV = 0,

дґ

^ + (У, У)у- +1 Ур' = 0, (1)

дґ р

•^2- + (У, У)р + рс2 divу' = 0,

дґ

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

^ + УЕ (О' ) = 0, (2)

о ґ

где О' = (, я2, Яз, я4) — вектор линеаризованных консервативных переменных,

Г (О') = (, Г) , Гр (О') р = х, у — векторы потоков, Гр (О') = Ар ()(3', О — вектор кон-

!~\ ЭГР (О)

сервативных переменных среднего течения, Ар --------матрица Якоби вектора потока

среднего течения. Заметим, что приведенное матричное уравнение (2) является более общей формой записи линеаризованных уравнений Эйлера, так как с формальной точки зрения включает возможность задания нестационарного и неоднородного фонового течения.

Следует отметить, что для гладких решений рассматриваемые далее методики могут быть обобщены и на случай полных уравнений Эйлера. При этом под конвективным потоком Г подразумевалась бы нелинейная функция полных консервативных переменных О = О + О'.

2. Повышение точности схемы по аналогии с высокоточными методами на структурированных сетках. Согласно методу контрольных объемов алгоритмы с определением переменных в узлах предусматривают построение дополнительных ячеек (контрольных объемов) вокруг узлов сетки. При этом в вычислительной газовой динамике (и аэроакустике в том числе) схемы конструируются путем аппроксимации законов сохранения на построенных элементарных ячейках.

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

| Е(О'(х, у, ґ))• пds. (3)

дСі

Здесь |Сг| — площадь расчетной ячейки, п — внешняя нормаль к границе ячейки ЭСг-.

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

Рис. 1. Барицентрический контрольный объем, или ячейка (слева), и двухсегментная грань контрольного объема, разделяющая точки 7 и к (справа)

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

Подобный подход можно применить и в случае нерегулярных сеток путем замены соответствующих центрально-разностных и направленных производных их «неструктурированными» аналогами. Полученная таким образом схема будет достигать своего высокого порядка, заданного одномерной аппроксимацией, только на структурированной сетке. При искажении структуры сетки точность алгоритма будет снижаться вплоть до второго, теоретически гарантированного, порядка в зависимости от качества сетки. Многопараметрическое семейство схем, обладающих вышеописанными свойствами, приводится в работах [1—4], а одномерный аналог данной схемы применительно к уравнению переноса подробно описан в [4]. Схемы третьего и пятого порядка точности из указанного многопараметрического семейства назовем МР83 и МР85 соответственно.

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

3. Повышение точности схемы за счет полиномиальной реконструкции переменных. Другим путем повышения точности схем на неструктурированных сетках является использование полиномиальной реконструкции переменных, заданных в узлах сетки, и уточнение на ее основе численного интегрирования в (3).

Пусть множество точек, соседних с точкой 7, есть множество О7. Тогда поток Г через границу ЭС7 ячейки С7 представим в виде:

| Г (О' (х, у, t ))• nds = ^ | Г (О' (х, у, t ))• nds =

дС7 ке°7 дС7 ПЭСк

2 (4)

= { Г(( ^ t))• n7kds,

ке°7 I=1 ЪС11к

2

где дС7 П дСк = ^ЭС7к — общая часть границы 7-й и к-й ячеек (см. рис. 1). Использование квад-

I=1

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

знания значений функции О' (х, у, ґ) в произвольных точках на границе ЭСг- П дСк. Но такая информация отсутствует, так как переменная О' = (, q'2, я3, q4 ) определена только в узловых точках. Эта проблема решается приближением исходных функций q1, q2, q3 , q4 некоторыми многочленами, определенными на ячейке С7, в том числе и на ее границе. Построение таких многочленов на примере абстрактной сеточной функции q(х, у, ґ) проводится следующим образом. Введем полином второго порядка Р" (х, у), приближающий сеточную скалярную функцию qn (х, у) в дискретный момент времени ґп на ячейке С7, так, чтобы выполнялось условие:

P

n (x, У) = q(x, У, tn), (5)

где РП (х, у) и яП — средние значения полинома и функции на ячейке, определенные как:

Р1П,1 (^ у) = итгIрп(x, у)dxdУ, я" = |7гт 1 Яп (x, у)dxdУ.

\С7\к \САк

Из равенства (5) следует вид полинома Ptn (x, y):

Pn (x, У T = I

p=1

5 ah

+ qi , (б)

где (X — Хо )р обозначает общий вид одночлена не выше второго порядка

( — Хо ) =( — Xо )(у — уо ), 0 <р + у< 2,

а (X — Х0 )а р — его среднее на ячейке С7

(X — Хо )р = С-1 (X — Хо )р dxdy. (7)

С7

Интеграл (7) для барицентрической ячейки вычисляется аналитически.

Далее определим среднее от полинома Р" ( у) по ячейке Ск (к Ф 7):

(S)

pn,k (x y )=Cr fpn (x y )dxdy=q+Z a, k, pap,,

Cklck p=l

A,k,p = C| Z f ( -X0 Tpdxdy -71 Z f ( -X0 )“

lCk| TeCk T \Ci\TeCiT

где T — треугольники, имеющие общую вершину в точках i или k соответственно.

По определению среднего полинома выполняется равенство P"k (x, y) = ql. Следовательно,

с учетом (8) имеем следующую систему линейных уравнений для определения неизвестных коэффициентов ap i:

Z A, k, pap, i = ql - qn. (9)

p=1

Рис. 2. Пример шаблонов N0: N0 =Оо (слева), N0 = Оо и и Ц (справа):

к=1

• — соседние первого уровня; Ф — соседние второго уровня

Пусть ячейки Ск являются соседними к ячейке С7 , т. е. к ёЦ . Если число соседних точек

Щ < 5, то система (8) недоопределена и имеет бесконечное множество решений. Если Щ > 5,

то система линейных уравнений переопределена и имеет единственное обобщенное решение в смысле минимальной невязки в среднеквадратичной норме (метод наименьших квадратов) . В связи с этим будем выбирать значения к из шаблона N., определенного как:

N =

ц , |Ц|

ц ц ии Ок к=1 |Ц|

Ц

где и Ц есть множество точек, соседних к точкам из множества Ц (рис. 2), иначе говоря, со-

к=1

седей второго уровня.

Таким образом, для определения неизвестных коэффициентов полинома (6) воспользуемся методом наименьших квадратов. Для этого потребуем минимального отклонения среднего значения данного полинома по всем соседним ячейкам, входящим в шаблон N., от средних значений функции, определенной в центрах** этих соседних ячеек, т. е. будем минимизировать целевую функцию:

кеN.

kеN7

Я7 — 2°р, Ак,р р=1

(10)

Из условия минимума целевой функции дН"/дап 7 = 0, я = 1,...,5 получаем систему линейных уравнений с квадратной невырожденной матрицей:

5

2 ар, ! Мр, п, 7 = К 7 , п = 1-5, р=1

kеN7

Мр, п, 7 2 ' ^7,к,р&7,к,п, |Мр, п, 7 Ц Ф 0, р, п . -5,

к., = 2 (

kеN7

Я 7 I), к, п.

: При Щ = 5 и йе^Д к р|| Ф 0 система (8) имеет единственное классическое решение. : Здесь под центрами ячеек подразумеваются узлы сетки, порождающие эти ячейки.

2

Решая систему (11) для каждой точки 7, получаем наборы значений коэффициентов полинома рп (x, у), аппроксимирующего значения сеточной функции яп (x, у). Заметим, матрицы М р, п, 7 не зависят от времени и могут быть обращены до начала выполнения численного интегрирования по времени.

Заменяя функцию О/п (x, у) = (п, , я'ъ , яТ ) ее полиномиальным приближением по опи-

санному выше методу, можно построить квадратуру Гаусса высокого порядка для вычисления поверхностного интеграла, входящего в равенство (4):

I Г((у)= I Г(п(^у))•) = 2®тГ(Рп(1,у1т))• )к. (12)

дС1, дС1, т=1

гк гк

Здесь ют и (т, у1т ), т = 1, ..., 5 — весовые коэффициенты и узловые точки квадратуры

Гаусса соответственно. Так как Рп (x,у) является полиномом второй степени, а квадратура Гаусса точна для полиномов степени 25 — 1, то имеет смысл выбрать значение 5, равное двум. При 5 = 2 имеем:

1

®1 =®2 = -,

4 = 1 (1—^)( — 4,)+x^, x' = 1^ +^ ](4—)+x^, (13)

у! = 2(1—^)( — у,)+у,, у2 = 1 (1 ^-13]( — у,)+,

где (, у^ | и , у,) — координаты вершин отрезка дС17к (см. рис. 2).

Таким образом, согласно (4), (12) и (13), поверхностный интеграл в выражении (3) вычисляется с третьим порядком точности. Следовательно, пространственная аппроксимация системы (2) также будет иметь порядок не ниже третьего вне зависимости от выбранного метода определения численного потока.

Зададим численный поток ф(кп, ОГ, п7к) через границу контрольного объема дСк с нормалью п7к как функцию двух переменных, удовлетворяющих условию согласования

ф(СОп, О п, п7к | = Г(п, п7к). Здесь Окп и О'п есть предельные значения функции О'п (x,у) на

границе дС, изнутри контрольных объемов С, и С7 соответственно. Так как функция О'п (x,у)

приближена полиномами в соответствующих ячейках, то численный поток определится по следующей формуле:

Г (Р( ^, у1т, t) )• )к =Ф(Р7 (т, у( , t), Рк (т, у1т, t), п!к ), (14)

где функция ф|^1, (2, п1, | может вычисляться, например, как поток с направленными разнос-

тями:

др(z1, ^2 )- п!к

(22 -^), 0<У< 1. (15)

+1, ] +1

'-1-7 ' ' 7

' I- / I '• / I

Рис. 3. Фрагмент «декартовой» сетки

Заметим, что именно такой способ определения численного потока использовался далее при проведении тестовых расчетов.

Выражения (3), (4), (12), (14), (15) определяют пространственную аппроксимацию системы уравнений (2), лежащую в основе схемы, обозначаемой в дальнейшем как БУ3.

Для проверки точности пространственной аппроксимации рассмотрим двумерное скалярное линейное уравнение переноса:

ди , ди

а— + Ь— = 0 дх ду

(16)

на треугольной декартовой сетке (рис. 3).

Применяя к уравнению (16) схему (3), (4), (12), (14), (15) и разлагая в ряд Тейлора, имеем следующие результаты: при у= 1

ди ди |

а-------+ Ь — = - а

дх ду

1 д4

12 дх4

и 3 II

Ах - а

1 д4и . 2 . , , 1

■|ь|

1 д4

12 дх4

и А 3

Ах

Ах Ау - а

д 4и

2 2 АхАу -2 дх ду 12 дх ду

1Ь1 ^АхАу2 - |Ь| дд2ди.2 Ах2Ау + 0 (■Ау4, Ау4) = 0 (Ау3, Ау3 ),

'12 дхду

'12 дх2ду2

при у= 0

ди ди

а-------+ Ь —

дх ду

= 0 (Ау4, Ау4).

Отсюда следует, что построенная схема действительно обладает порядком точности пространственной аппроксимации не ниже третьего. В случае центрально-разностного аналога схемы (у= 0) порядок аппроксимации — четвертый.

Перед тем, как перейти к обсуждению тестовых расчетов, отметим, что в схемах обоих рассмотренных классов интегрирование по времени для полудискретной системы (3) может осуществляться одним из методов повышенной точности, например, по «классической» схеме Рунге — Кутты четвертого порядка, которая и была использована в настоящей работе.

4. Примеры тестовых расчетов. В качестве теста, на котором на практике изучались свойства предложенных схем, рассматривалась начальная задача о распространении гауссового импульса, заданного возмущениями плотности и давления, в неподвижной среде:

1д2/ 2

р(х, у,0) = р(х, у, 0) = е Ь

(х2 + у2)

' (х, у,0) = V (х, у,0) = 0,

где Ь — полуширина гауссового импульса. Заметим, что такая постановка задачи обеспечивает наличие в решении только акустических волн.

Задача решалась в области 100 X 100 на двух сетках с характерными размерами ячеек 1

и у Фрагменты трех типов сетки, используемых в расчетах, изображены на рис. 4. Численное решение сравнивалось с точным при значении времени 40, порядок точности в нормах

Ьр, р = 1,2 оценивался по возмущениям плотности согласно формуле:

Ит

1о§>|\рк=1 рехас1;|1°§>| |р/г=12 рехас1;|

1о§ (2 )

Рис. 4. Фрагменты сеток, участвующих в тестовых расчетах (слева направо): «декартовая», «остроугольная», «тупоугольная»

Задачей тестирования было сравнение схем МР83, МР85 третьего и пятого порядка в точности многопараметрического семейства (см. раздел 2) и схемы БУ3 третьего порядка точности, полученной с использованием полиномиальной реконструкции переменных (см. раздел 3). Напомним, что под точностью схем МР83 и МР85 имеется в виду точность аппроксимации, теоретически достигаемая на треугольной декартовой сетке.

В качестве примера изолинии возмущений плотности, полученные при помощи схемы БУ3 на «остроугольной» сетке (см. рис. 4) размерностью 39 527 узлов для случая полуширины гауссового импульса Ь= 6, приведены на рис. 5. Рис. 5 показывает 10 изолиний в диапазоне от -0.0324 до 0.0580 с равномерным шагом. Полные данные по оценкам поряд-

ка точности Ыь схем для различной полуширины

Рис. 5. Задача о распространении гауссового импульса в неподвижной среде при начальной полуширине импульса Ь = 6 (изолинии возмущений плотности на момент времени 40, полученные при помощи схемы ЕУ3)

гауссового импульса представлены в таблице.

Собранные в таблице данные показывают характерные, в том числе и несколько неожиданные,

свойства исследуемых схем. Следует сразу отметить, что существенно более низкие порядки для расчетов с Ь= 3 объясняются специально выбранным плохим разрешением по сетке начального возмущения (для грубой сетки с шагом 1 — три точки на полуширину гауссового импульса).

Полуширина импульса Ь = 3 Полуширина импульса Ь = 6

Схема ЕУ3 МРБ3 МРБ5 ЕУ3 МРБ3 МРБ5

Норма Сетка і («декартовая»)

іі 1.42 1.43 2.80 2.45 2.50 4.48

і2 1.23 1.25 2.45 2.32 2.38 4.35

Норма Сетка 2 («остроугольная»)

іі 1.57 1.55 3.07 2.68 2.70 2.98

і2 1.41 1.4 2.87 2.60 2.62 2.89

Норма Сетка 3 («тупоугольная»)

іі 1.33 1.22 1.98 2.66 2.44 2.04

і2 1.16 1.10 1.78 2.43 2.33 1.94

В то же время, эти результаты показывают явное преимущество в точности схемы MPS5 даже в случае «недекартовых» сеток, где они имеют теоретически только второй порядок. При лучшем разрешении начального импульса (b= 6) для всех используемых сеток схема FV3 подтверждает свой близкий к теоретическому, третьему, порядок аппроксимации. И в этом случае для сеток, отличных от «декартовых», схемы MPS3 и MPS5 показывают точность, заметно превышающую теоретическую, т. е. второго порядка. Исключением из наблюдаемой тенденции являются результаты для схемы MPS5, полученные на сетке, состоящей только из тупоугольных треугольников. Вероятным объяснением этого факта является «плохое» для схемы MPS5 качество сетки, определяющее «неудачный» для нее расширенный шаблон аппроксимации.

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

Работа выполнена при поддержке Российским фондом фундаментальных исследований (проекты № 09-01-00833-а, 09-01-12022-офи_м).

ЛИТЕРАТУРА

1. Debiez C. and Dervieux A. Mixed element volume MUSCL methods with weak viscosity for steady and unsteady flow calculation // Computer and Fluids. 1999. V. 29, pp. 89—118.

2. Abalakinl. V., Dervieux A. and Kozubskaya T. K. A vertex centred high order MUSCL scheme applying to linearized Euler acoustics. — INRIA report RR4459, April 2002, http://www. inria. fr/rrrt/rr-4459.html.

3. Abalakin I. V., Kozubskaya T. K., Dervieux A. High accuracy finite volume method for solving nonlinear aeroacoustics problems on unstructured meshes // Chinese J. of Aeroacoustics. 2006. V. 19, N 2.

4. Абалакин И. В., Козубская Т. К. Многопараметрическое семейство схем повышенной точности для линейного уравнения переноса // Математическое моделирование.

2007. Т. 19, № 7, с. 56—66.

5. Zhang Y.-T., Shu C.-W. Third order WENO scheme on three dimensional tetrahedral meshes // Commun. Comput. Phys. 2009. V. 5, pp. 836—848.

Рукопись поступила 18/III2009 г.