CC BY
82
Том XLI
УЧЕНЫЕ ЗАПИСКИ ЦАГИ 2010
№ 2
УДК 519.63:517.96 534.2
О РЕАЛИЗАЦИИ ЧИСЛЕННОГО МЕТОДА, СОХРАНЯЮЩЕГО ДИСПЕРСИОННОЕ СООТНОШЕНИЕ,
НА ДИХОТОМИЧЕСКИХ СЕТКАХ
В. Г. БОБКОВ, Д. С. ЗВЕНКОВ
Рассмотрена реализация DRP схемы, сохраняющей дисперсионное соотношение, на локально измельченных сетках и, в частности, различным типам интерполяции, сохраняющим свойства схемы. Рассмотрены оригинальные методы построения локально сгущающихся сеток, методы интерполяции решения на границах блоков с разным шагом сетки, а также обработка граничных условий.
Ключевые слова: схема, сохраняющая дисперсионное соотношение, DRP-интерполяция, искусственные граничные условия, локально измельченная сетка, дихотомическая сетка.
Расчет с высоким порядком точности в задачах вычислительной гидродинамики, а особенно в задачах вычислительной аэроакустики, требует достаточно точного разрешения различных пространственных масштабов. Качество решения вычислительной задачи в большой степени зависит от качества моделирования областей, в которых течение имеет особенности и, в частности, включает мелкомасштабные структуры, влияющие на течение в целом. Для разрешения этих особенностей необходимо наличие достаточно подробной сетки в упомянутых выше областях.
Одним из подходов к решению этой проблемы является использование расчетных сеток, локально сгущенных в областях, требующих особо высокого разрешения. Так, использование сеток, сгущенных специальным «дихотомическим» образом (здесь и далее под «дихотомической» сеткой имеется в виду блочная декартовая локально сгущенная сетка, в которой линейные размеры соседних ячеек, если они разные, могут различаться только в 2 раза) совместно с численными методами повышенной точности, позволяет получить положительные результаты в задачах вычислительной аэроакустики за счет:
локально-декартового характера блоков сетки одного уровня, позволяющего применять метод высокого порядка точности для декартовых сеток;
оптимальной структуры хранения сетки на основе четверичных/восьмеричных деревьев (quadtree/octree);
возможности адаптации сетки для разрешения разномасштабных структур; естественного распространения подхода на мультипроцессорные платформы.
Подобные сетки не являются уникальными, они упоминаются в исследовательских работах [1] и используются в таких вычислительных пакетах, как FlowVision [2] и NUMECA [3]. Заметим, что в упомянутых работах на подобных сетках построены численные схемы второго порядка, в то время как в данной работе представлен алгоритм четвертого порядка аппроксимации по пространству.
Таким образом, помимо оригинального метода построения дихотомической сетки здесь рассмотрен пример реализации метода DRP (Dispersion Relation Preserving) [4], сохраняющего дисперсионное соотношение, на сетках описанного выше типа применительно к задачам вычислительной аэроакустики. В частности, затронуты проблемы интерполяции решения на границах блоков с разным шагом сетки, локального интегрирования по времени и обработки граничных условий на примере методик, реализованных в комплексе программ WHISPAR [5, 6].
Построение сетки. Как уже было упомянуто выше, в данной работе рассматривается прямоугольная локально-измельченная сетка, построенная с использованием четверичного (или восьмеричного — в трехмерном случае) дерева. В таком представлении каждой расчетной ячейке сетки соответствует лист дерева, а совокупность всех листьев дерева полностью определяет всю расчетную сетку.
Такое представление сетки дает несколько ключевых преимуществ перед другими возможными вариантами обработки и хранения сетки. Одним из важных преимуществ является быстрая и эффективная процедура измельчения и загрубления сетки, которая может применяться как при построении сетки для стационарных задач, так и для динамического изменения сетки в процессе расчета нестационарных задач.
Построение сетки происходит путем последовательного измельчения («роста») основного дерева, т. е. путем разбиения выбранных ячеек сетки, что эквивалентно порождению новых потомков у соответствующих листьев исходного дерева (рис. 1).
При этом, с каждым листом дерева (т. е. с каждой ячейкой сетки) ассоциирован его уровень — число делений, с помощью которых была получена данная ячейка.
Как правило, измельчение сетки необходимо в пристеночных областях, где существенное влияние на общую картину течения оказывает точность моделирования газодинамических эффектов вблизи стенки.
Основная идея рассматриваемого метода построения сетки заключается в построении сгущения сетки в соответствии с наперед заданной функцией / (г )е [0; 1], г еО плотности сетки.
Эта функция каждой точке расчетной области ставит в соответствие желаемый уровень измельчения сетки следующим образом: нулевое значение функции соответствует минимальному («фоновому») уровню измельчения сетки, а 1 — максимальному.
Этот алгоритм может быть дополнен функцией-маркером, определяющей реальную геометрию твердых тел От в расчетной области:
Л(г ) =
[0, г е О/От, [ 1, г е От.
Собственно построение сетки с использованием описанной методики выглядит следующим образом (пример построения сетки вокруг цилиндра — см. рис. 2):
1) изначально создается так называемая «фоновая» сетка — равномерное заполнение расчетной области ячейками минимального требуемого уровня;
Рис. 1. Пример локально-измельченной сетки и соответствующее ей четвертичное дерево
X
Рис. 2. Пример построения локально-сгущенной сетки вокруг цилиндра
2) с использованием функции плотности сетки для каждой ячейки производится проверка, не достигнут ли максимальный уровень для данной ячейки;
3) если максимальный уровень не достигнут, производится деление ячейки и процедура с шага 2 рекурсивно вызывается для новых полученных ячеек;
4) с использованием функции-маркера /т (г) расчетные ячейки, лежащие внутри твердых
тел, помечаются как «мертвые».
Такой подход к построению сетки позволяет строить сетки для достаточно сложных конфигураций и имеет ряд весьма полезных свойств:
1) построение сетки сводится к построению функций, описывающих плотность сетки (f (r))
и геометрию твердых тел (т (г)), не зависящих от текущих свойств сетки;
2) в случае наличия в исходной задаче нескольких требований к измельчению сетки и множества твердых тел очевидным образом могут быть построены и использованы суперпозиции соответствующих функций:
f (г ) = max[ fi (г^..^ fK (г)], FеЦ
/т () = max[ fT l (),..., fT,L ()], г е Q;
3) функция-маркер может быть использована для уточнения исходной функции плотности сетки для дополнительного измельчения сетки вблизи твердых границ;
4) все описанные алгоритмы применимы для решения как двумерных, так и трехмерных задач.
Предложенная методика дает гибкий и эффективный инструмент для построения дихотомических сеток, наиболее точно описывающих геометрию задачи с заданным уровнем детализации сетки.
DRP-интерполяция. В данной работе в качестве основной численной схемы для вычисления пространственных производных использовалась схема, сохраняющая дисперсионное соотношение (Dispersion Relation Preserving), предложенная К. Тамом [4]. Выбранная схема обладает достаточно широким симметричным пространственным шаблоном, во многом определяющим специфику предложенных методик и алгоритмов.
В основе используемой схемы лежит аппроксимация производной функции на семиточечном пространственном шаблоне:
Ї(х “>-г (х+1 Ах >•
1 3
где коэффициенты а. выбираются на основании двух условий: условия четвертого порядка
аппроксимации (путем разложения обеих частей уравнения в ряд Тейлора и выбора условий аппроксимации для коэффициентов) и условия минимизации разницы между образами Фурье левой и правой частей:
іа/--
( і 3 — У а,вЩА
Ах ^3 1
V 1=-3
Л п/ 2
/, Е = Г |аАх-аАх\ё(аАх)
-п/ 2
>Ш1П,
а1е
і а] Ах
— эффективное волновое число. Полученная таким образом схема поми-
где а = — У Ах
]=-3
мо четвертого порядка аппроксимации по пространству обладает улучшенными дисперсионными характеристиками, что позволяет с успехом использовать ее в расчетах аэроакустических задач. Однако ограничение на использование с данной схемой исключительно центральных симметричных шаблонов существенно снижало ее гибкость как вычислительного инструмента и возможность применения ее на описанных выше сетках. В связи с этим было решено несколько модифицировать используемые вычислительные алгоритмы.
□ Е П ■ В С ] □ А ■ п ■ 1 1 1 1 1
■ 1 1 1 1 1
ЫНЫ т т т
Рис. 3. Пространственный шаблон вблизи интерфейса между блоками разного уровня (в точках А, В, С необходимо интерполировать значение функции)
Так для расчета на описанных выше дихотомических сетках были разработаны методики БКР-интерполяции [7, 8], используемой для расчета численных производных вблизи интерфейсов блоков сетки с разными шагами по пространству, где вычислительный шаблон схемы, построенный на существующих точках сетки, является несимметричным (рис. 3). Вблизи интерфейса блока с большим шагом имеются точки (А, В, С на рис. 3), в которых тем или иным способом необходимо восстановить значение функции, чтобы в последующем вычислить производную, используя симметричный шаблон. Интерполировать значения в этих точках можно различными способами, причем крайне важно качество интерполяции. Авторами был исследован ряд алгоритмов, начиная с билинейной и заканчивая бикубической и сплайн-интерполяцией, но качество результата было неудовлетворительным — в процессе численного решения на интерфейсах блоков возникали нефизичные высокочастотные осцилляции. Наиболее качественный результат дала интерполяция, сконструированная на основе БКР-подхода: интерполируемое значение функции вычислялось с использованием шеститочечного симметричного шаблона по следующей формуле:
/(х)« /(х) = У (х + (. - 2.5)Ах)
у=о
где коэффициенты Sj получены в результате минимизации функционала
—K
i (j—2.5) aAx
Sje J =1
■У
d aAx,
который представляет собой разницу между точным и интерполированным значениями для Фурье гармоник из определенного диапазона [—K, K] волновых чисел.
Заметим, что при вычислении значений в интерполяционных точках могут использоваться вновь интерполированные значения, в результате чего при вычислении производной задействуется до 36 точек. Реализация алгоритма представляет собой рекурсивную процедуру. Так, если для вычисления интерполированного значения функции необходима точка, не принадлежащая исходной сетке, рекурсивно вызывается эта же процедура интерполяции для недостающей точки и т. д. до тех пор, пока интерполяция не будет использовать только точки заданной сетки.
Параллельная реализация. Естественным путем распараллеливания для явного DRP-метода является использование геометрического параллелизма, основанного на разбиении сетки на части. При этом каждая подобласть сетки обсчитывается на своем вычислительном узле. Рассматриваемая DRP-методика вычисления производных на интерфейсах блоков сетки вызывает сложности в реализации распределенной версии. Трудность вызвана тем обстоятельством, что до вызова рекурсивной процедуры вычисления интерполированного значения неизвестно, какие точки необходимы для вычисления производной в данной точке, и, следовательно, нет возможности определить, данные каких точек необходимо пересылать соседним вычислительным узлам. Однако эта проблема имеет решение — достаточно на заданной сетке один раз взять производные по всем направлениям и в процессе интерполирования пометить точки, задействованные при численном дифференцировании. Таким образом, все данные, необходимые для определения зон пересылки между вычислительными узлами, станут известны. В то же время, недостатком такого алгоритма, снижающим эффективность его распараллеливания, является заметное различие от точки к точке в объемах пересылаемой информации. С этой точки зрения наилучший результат могут давать способы разбиения сетки на подобласти, граница между которыми минимальным образом пересекается с границами между блоками сетки с разными пространственными шагами. Возможно, решить эту проблему удастся с помощью альтернативных методов упорядочивания ячеек и последующего разбиения сетки на основе кривых, заполняющих пространство. Авторами была предпринята попытка использования разбиения на основе кривой Гильберта [9]. Однако, несмотря на ряд положительных свойств (таких, как крайне компактное представление сетки, быстрый поиск ячеек сетки и близкое к оптимальному (в смысле количества ячеек в каждой подобласти) разбиение), качество разбиения оказалось хуже, чем при использовании классических алгоритмов разбиения сетки на основе графов связей (библиотека METIS). Это связано с тем, что методика на основе кривой Гильберта наилучшим образом подходит для конечно-объемного подхода, в то время как в данной работе используются конечно-разностные схемы. Поэтому в настоящее время авторами ведется разработка конечно-объемной модификации схемы DRP для последующей апробации предложенного метода разбиения сетки с использованием кривой Гильберта.
Результаты. В качестве тестовой проблемы была выбрана задача о рассеянии точечного возмущения на наклонной стенке [10] в несколько усложненной постановке. В расчетной области (—100; 100) х(0; 200) в точке (0; 40) в начальный момент времени задавалось возмущение плотности и давления в виде гауссова импульса следующего вида:
—ln 2
Р(, y, t = 0) = p(, y, t = 0) = e
x 2+(y—25)2 25
Нижняя часть расчетной области была ограничена наклонной твердой стенкой, проходящей через нижнюю правую точку области (100; 0) и через точку (-100; 20) (рис. 4). Данная задача
имеет аналитическое решение [10].
Для численного решения задачи использовались лианеризованные уравнения Эйлера. На поверхности стенки ставились хорошо себя зарекомендовавшие отражающие граничные условия Тама [11], на свободных границах — выходные неотражающие граничные условия Тама [12].
X
Рис. 4. Сетка в тестовом расчете (крупная)
X
Рис. 5. Результат расчета тестовой задачи на мелкой сетке (изолинии
плотности)
Пример расчета приведен на рис. 5, на котором изображены изолинии плотности в момент
времени t = 50. В таблице для поля плотности приведены отклонения численного решения от
точного в норме L2 на разные моменты времени, а также фактический порядок схемы, вычис-
, и икрупн , и ммелк
l0g| |p—pт0ч^|| L ~ ^0g| |p—PTOчHI L2
ленный по формуле П =-------------;----2--;------;------;-----.
log ((крупн) — log (мелк)
Время, t L2 норма, сетка 2h L2 норма, сетка h П
10 5.098E-001 3.879E-002 3.72
20 5.304E-001 4.211E-002 3.65
30 5.462E-001 4.410E-002 3.63
40 5.588E-001 5.364E-002 3.38
50 5.762E-001 5.479E-002 3.39
В таблице приведены разницы между точным решением и численным, полученным на сетке, сгущающейся к наклонной стенке и имеющей три уровня сгущения, причем данные приведены для двух сеток — крупной (Ahmin = 1.6) (см. рис. 4.) и мелкой (Ahmin = 0.8).
Заключение. В процессе работы были изучены, реализованы и апробированы перспективные численные методики. Среди них: оригинальный алгоритм построения качественной локально измельченной сетки дихотомного типа, метод высокоточного вычисления DRP-производных на вышеупомянутых сетках с использованием DRP-интерполяции на границах интерфейсов блоков сетки с разным пространственным шагом, подход к моделированию тел с твердыми границами криволинейной конфигурации путем автоматического сгущения сетки вблизи границ и др.
Описанные методики и алгоритмы были реализованы в комплексе программ WHISPAR, что позволило проводить высокоточные расчеты задач аэроакустики в условиях достаточно сложной геометрии. Дальнейшая работа будет направлена на улучшение реализации свободных граничных условий, верификации комплекса программ на серии тестовых задач и адаптации реализованных методик и алгоритмов для трехмерного случая.
ЛИТЕРАТУРА
1. ChetverushkinB., ChurbanovaN., Trapeznikova М., Sukhinov A., Malinovskij A. Adaptive cartesian mesh refinement for simulating multiphase flows in porous media // Computational Methods in Applied Mathematics. 2008. V. 8, N 2, p. 101 —115.
2. Аксенов А. А., Гудзовский А. В. Пакет прикладных программ Flow Vision // Труды МФТИ, сер.: Аэрофизика и прикладная математика. 1998, с. 45—56. URL: http://www.flowvision.ru (дата обращения 02.02.2009).
3. Официальный сайт вычислительного пакета NUMECA. URL: http://www.numeca.com (дата обращения 02.02.2009).
4. Tam C. K. W., Webb J. C. Dispersion-Relation-Preserving finite difference schemes for computational acoustics // J. of Computational Physics. 1993. V. 107, N 2, p. 262—281.
5. Звенков Д. С. О реализации схемы DRP на блочных сетках // Математическое моделирование. 2007. Т. 19, № 7, с. 112—118.
6. БобковВ. Г. Краткое описание программного комплекса WHISPAR для моделирования задач аэроакустики на структурированных сетках // Математическое моделирование.
2007. Т. 19, № 8, с. 123 — 128.
7. Tam C. K. W., Kurbatskii K. A. Multi-Size-Mesh Multi-Time-Step Dispersion-Relation-Preserving scheme for multiple-scales aeroacoustics problems // J. of Computational Fluid Dynamics. 2003. V. 17, N 2, p. 119—132.
8. Tam C. K. W., Kurbatskii K. A. A wavenumber based extrapolation and interpolation method for use in conjunction with high-order finite difference schemes // J. of Computational Fluid Dynamics. 2000. V. 157, N 2, p. 588—617.
9. Bobkov V., Kozubskaya T., Zvenkov D., Tsikaridze Sh. On locally refinement mesh processing for parallel CFD // Proceedings of International Conference NUMGRID-2008. 2008, p. 93—97.
10. Hardin J. C., Ri store11i J. R., T am C. K. W. ICASE/LaRC workshop on benchmark problems in computational aeroacoustics (CAA) // Hampton: NASA Conference Publication-3300, 1994, 438 p.
11. Tam C. K. W., Dong Z. Wall boundary conditions for high-order finite difference schemes in computational aeroacoustics // Theoretical and Computational Fluid Dynamics. 1994.
V. 6, N 5—6, p. 303—322.
12. Bayliss A., Turkel E. Radiation boundary conditions for wave-like equations // Communications on Pure Applied Mathematics. 1980. V. 33, N 6, p. 707—725.
Рукопись поступила 25/III2009 г.
