Журнал вычислительной математики и математической физики, 2022, T. 62, № 10, стр. 1620-1631

ВВЕДЕНИЕ

Задача обтекания подвижных тел имеет исключительно важное значение для современной вычислительной аэродинамики. Интерес в одновременном рассмотрении двух или более обтекаемых объектов, каждый из которых движется по своим законам, стимулирует развитие различных методов моделирования и техники соответствующей работы с геометрией. В данной работе предлагается новая методика сеточной анизотропной адаптации с целью ее дальнейшего использования в комбинации с методом погруженных границ при решении задач обтекания тел сложной формы. Среди алгоритмов погруженных границ в настоящей статье мы ориентируемся на метод штрафных функций Бринкмана [1], [2], в рамках которого обтекаемый объект моделируется внесением дополнительных источниковых членов в газодинамические уравнения, в качестве которых мы подразумеваем уравнения Эйлера, Навье–Стокса и осредненные по Рейнольдсу уравнения Навье–Стокса для описания течений вязкого сжимаемого газа. Источники в уравнениях определяют тело как сплошную среду со слабой проницаемостью, их величины отличны от нуля только в точках, лежащих внутри геометрического положения движущегося объекта. Таким образом, метод погруженных границ позволяет определять задачу внешнего обтекания в односвязной области, что, в свою очередь, открывает возможность использования техники перераспределения узлов при сохранении сеточной топологии для уточнения локализации поверхности обтекаемых объектов.

Для расчета задач внешнего обтекания при “классическом” подходе, т.е. на сетках, граничные узлы (или элементы) которой согласованы с границей тела, существует ряд известных методов. Среди них сеточная деформация [3] является наименее вычислительно затратным алгоритмом. Она позволяет достаточно быстро менять геометрию сетки для нового положения тела на каждом шаге по времени. Недостатком сеточной деформации является тот факт, что при больших амплитудах движения качество сетки становится слабо контролируемым вплоть до разрывов и перепутываний, а потому для больших смещений этот подход не годится. Альтернативный подход, основанный на перестроении сетки [4], как локальном, так и глобальном, в том числе и на доразбиении [5], хотя и может быть достаточно быстрым, однако требует вовлечения в численный алгоритм нежелательной интерполяции, а также заметно усложняет процесс распараллеливания, так как топология и размер сетки могут со временем меняться. Метод Химеры [6], подразумевающий наличие двух или нескольких наложенных сеток, также требует интерполяции и усложненной реализации обмена расчетных данных между сетками.

В статьях [7]–[10] разрабатываются и исследуются алгоритмы для адаптации сетки с использованием управляющей функции или метрик. Существуют работы по подбору метрики и для анизотропных сеток как в контексте генерации или доразбиения сетки [5], так и для подвижных адаптивных сеток [11]. Способы задания различных управляющих метрик в различных методах адаптации обсуждаются в работе [12]. Как правило, все представленные в литературе алгоритмы подразумевают оценку апостериорной ошибки, на основе которой и строится метрический тензор.

Следует отметить, что адаптация подвижной сетки в задачах вычислительной газовой динамики при использовании метода погруженных границ использовалась и ранее (см., например, работы [13], [14]). Отличительной чертой предлагаемой нами методики сеточной адаптации является, как уже отмечено выше, возможность использования техники перераспределения узлов с сохранением топологии исходной сетки. Это свойство адаптации позволяет избежать снижения точности численного результата из-за интерполяции, сохраняет структуру алгоритма и, как следствие, сильно упрощает реализацию параллельной версии. Одной из трудностей при таком подходе применительно к реальным промышленно-ориентированным задачам для высокоскоростных течений является обеспечение необходимого сеточного разрешения в пристеночной области. При этом требуемое разрешение должно сохраняться независимо от возможных перемещений обтекаемого тела. При построении подхода мы предполагаем, что исходная сетка весьма умеренных размеров в начале расчета сгущается к поверхности “погруженного” в расчетную область объекта, а затем начинает передвигаться вместе с ним. Желаемая степень сеточного сгущения к поверхности тела при этом определяется требованием размера шага сетки в нормальном направлении внутри пристеночного пограничного слоя. Толщина же пограничного слоя уменьшается с ростом числа Рейнольдса, характеризующим уровень турбулентности течения. Это приводит к жесткому ограничению на малость шага сетки вблизи поверхности. В настоящей работе, посвященной только методу адаптации как таковому, не исследуется вопрос о достижимости нужного разрешения сетки в пограничном слое. Безусловно, ответ на этот вопрос зависит от количества узлов исходной сетки и, в общем случае, от их начального распределения, т.е. от тех параметров, которые можно регулировать на стадии препроцессинга. При этом стоит заметить, что, в случае недостатка узлов и недостаточного разрешения погранслойной области сеткой, при моделировании турбулентных течений мы предполагаем использовать технику пристеночных функций.

Практически во всех методах адаптации к поверхности тела, в том числе и в нашем подходе, главным входным параметром для адаптации выступает функция расстояния. Понятно, что в задачах с подвижными объектами скорость вычисления расстояния в узлах динамически изменяемой сетки становится критической с точки зрения обеспечения эффективности всего алгоритма в целом. Для минимизации вычислительной стоимости этой процедуры мы предлагаем хранить все параметры геометрии обтекаемого тела и, в первую очередь, поле расстояния до его поверхности в узлах восьмеричного дерева, называемом далее интерполяционной решеткой [15]. Такое дерево является атрибутом тела и строится один раз на этапе подготовки геометрии. Далее интерполяционная решетка движется вместе с телом и позволяет вычислить новое расстояние в любой точке расчетной области за конечное время. Такая методика уже тестировалась на различных конфигурациях и использовалась при численном решении двумерных задач внешнего обтекания подвижных тел [16], [17]. Решение же существенно трехмерных задач вычислительной газовой динамики потребовало дальнейшего развития предлагаемого подхода. Дело в том, что изотропный характер ранее предложенной адаптации, определяемой лишь расстоянием до поверхности тела, не дает возможности построения сильно анизотропных сеток вдоль поверхности обтекаемого тела, характерных для численного решения широкого класса задач внешнего обтекания с высокими числами Рейнольдса. Действительно, изотропное сгущение сетки к поверхности, обеспечивающее требуемую малость пристеночного шага, с неизбежностью приводит к построению сеток чрезмерно большой размерности даже в двумерном случае и неприемлемо большой размерности – в трехмерном. Это обстоятельство потребовало учета в разработанном алгоритме анизотропии геометрии. Для качественной реализации анизотропной адаптации для тел сложной формы, особенно в трехмерной постановке, помимо функции расстояния, необходимо учитывать информацию более высоких порядков для входящей геометрии. В настоящей работе описывается наша реализация метода учeта кривизны при построении управляющей метрики адаптации. Кривизна рассчитывается для триангулированной поверхности подвижного тела, а ее значения записываются в вершины интерполяционной решетки. В процессе сеточной адаптации эта информация используется для настройки управляющих параметров метрики. В статье методика анизотропной адаптации формулируется как для двумерной, так и трехмерной постановки.

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

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

1. МЕТОД АДАПТАЦИИ ПОДВИЖНОЙ НЕСТРУКТУРИРОВАННОЙ СЕТКИ

В методе адаптации подвижной сетки расчетная сетка рассматривается как упругий материал, который подвергается сжатию вблизи границы движущегося тела. Упругая деформация, зависящая от времени, задается d-мерным отображением $x(\xi ,t)\,:{{R}^{d}} \times R \to {{R}^{d}} \times R$.

Пусть $C$ обозначает матрицу Якоби отображения $x(\xi ,.)$, где ${{c}_{{ij}}} = {{\partial {{x}_{i}}} \mathord{\left/ {\vphantom {{\partial {{x}_{i}}} {\partial {{\xi }_{j}}}}} \right. \kern-0em} {\partial {{\xi }_{j}}}}$. Пусть ${{x}^{n}}(\xi )$ – отображение исходной сетки в сетку на момент времени ${{t}^{n}}$.

Чтобы сформулировать вариационную задачу для сеточной адаптации на дифференциальном уровне, нужно ввести предположение, что ${{x}^{n}}(\xi )$ – квазиизометрический диффеоморфизм. Деформация на момент времени ${{t}^{{n + 1}}}$ определяется решением вариационной задачи [18]

В двумерном случае мы рассматриваем отображение треугольника с координатами вершин $\left( {{{h}_{0}},{{h}_{1}},{{h}_{2}}} \right)$ в лагранжевой системе координат в области ${{\Omega }_{\xi }}$ на треугольник с координатами вершин $\left( {{{p}_{0}},{{p}_{1}},{{p}_{2}}} \right)$ в эйлеровой системе координат в области ${{\Omega }_{x}}$. Лагранжевы координаты вморожены в начальную сетку, а эйлеровы являются декартовыми координатами в расчeтной области. В данном случае области $\Omega = {{\Omega }_{x}} = {{\Omega }_{\xi }}$. Тогда Матрица ${{G}_{x}}(x,t) = {{Q}^{{\text{т}}}}Q$ определяет метрический тензор в эйлеровых координатах. В функционале (1) функция $W(C)$ обозначает поливыпуклый упругий потенциал (внутреннюю энергию), который является взвешенной суммой меры искажения формы и меры искажения объема: В большинстве рассматриваемых случаев использовалось значение параметра $\theta = 4{\text{/}}5$. Упругий потенциал достигает минимума, когда деформация изометрична, т.е., когда допустимы только поворот и параллельный перенос.

Пусть граница подвижного деформируемого тела в момент времени t задается в виде нулевой изоповерхности скалярной функции $u(x,t)$. Предполагается, что $u(x,t)$ близка к функции расстояния со знаком от границы тела, поскольку в реальных задачах вычисление точной функции расстояния, как правило, не имеет смысла, а используются ее различные приближенные варианты.

Метрический тензор ${{G}_{x}}$ может быть записан в виде ${{G}_{x}} = U{{\Sigma }^{2}}{{U}^{{\text{т}}}}$, где столбцы $U$ являются собственными векторами ${{G}_{x}}$, а ${{\Sigma }^{2}}$ – диагональная матрица с элементами на диагонали $\sigma _{i}^{2}$, т.е. собственными значениями ${{G}_{x}}$. В каждой точке расчетной области можно задать локально оптимальные координаты с помощью аффинного преобразования с матрицей $Q = \Sigma {{U}^{{\text{т}}}}$. Говорят, что расчетная сетка оптимальна (изометрична), если после, применяя данное аффинное преобразование к ячейкам полученной сетки, получится начальная сетка, или

где V – произвольная ортогональная матрица.

Метрика адаптации $G = {{G}_{x}}$ требует задания нормального и тангенциального направлений в каждой точке $p$ (столбцы матрицы $U$), которые в свою очередь определяются изоповерхностью функции $u(x,t)$, проходящей через $p$ , а также коэффициентов растяжения ${{\sigma }_{i}}$ вдоль этих направлений. Чем больше значение ${{\sigma }_{i}}$, тем меньше размер ячейки в направлении столбца ${{u}_{i}}$. Соотношение ${{\max }_{{i,j}}}\frac{{{{\sigma }_{i}}}}{{{{\sigma }_{j}}}} - 1$ определяет степень анизотропии. Функция ${{\sigma }_{1}} = {{\sigma }_{{{\text{normal}}}}}(x,t)$ определяет растяжение сетки вдоль направления, нормального к телу, а ${{\sigma }_{2}} = {{\sigma }_{{{\text{tangential}}}}}(x,t)$ (${{\sigma }_{{2,3}}}$ в 3D) определяет пространственное распределение анизотропии. В узком слое на границе тела анизотропия максимальна и стремится к нулю вдали от тела.

Очень часто некорректно построенная метрика может приводить к резкому росту константы квазиизометрии, скачкам размеров ячеек и визуальным “разрывам” в сетке. Оптимальная метрика, которая дает наименьшую константу квазиизометрии, не известна даже для простых тел. В данной работе используются некоторые эвристические правила задания метрики, когда растяжение ${{\sigma }_{1}}$ принимается константой в слое на границе тела и убывает по гиперболе при отдалении от него, а тангенциальное растяжение определяется кривизной поверхности и постепенно уменьшает свое влияние вдали от границы. Для тела задается некоторое пятно радиуса $D$, за пределами которого метрика ${{G}_{x}}$ равна единичной матрице.

Метрика изотропной адаптации определяется соотношением

где ${{\sigma }_{1}}(x,t) = \phi \left( {u(x,t)} \right)$ – сеточная плотность. В данном случае степень анизотропии равна нулю. Использование изотропной метрики во всей расчетной области возможно только для небольших значений ${{\sigma }_{1}}$, так как решением вариационной задачи в этом случае является отображение, при котором прообраз тела почти в ${{\sigma }_{1}}$ раз больше размера тела в эйлеровой системе координат. Это приводит к стягиванию внутрь тела большого числа точек расчетной области.

Анизотропная версия метрического тензора определяется соотношением

На участках большой кривизны границы или около острых углов задается соотношение ${{\sigma }_{2}} = {{\sigma }_{1}}$, в остальной части границы задано соотношение ${{\sigma }_{2}} = {{\sigma }_{1}}{\text{/}}K$, где $K > 1$ – степень анизотропии, заданная извне. При отдалении от тела ${{\sigma }_{2}}$ приближается к ${{\sigma }_{1}}$, а дальше они оба стремятся к единице. Более подробное описание задания ${{\sigma }_{1}}$ и ${{\sigma }_{2}}$будет приведено далее.

К функционалу (1) применяется стандартная конечно-элементная дискретизация. Метрический тензор ${{G}_{{}}}$ хранится в узлах сетки, что критично для стабильности тонких и сильно сжатых сеточных слоев. Для решения задачи оптимизации на каждом шаге по времени применяется техника градиентного спуска с предобуславливателем [19], где направление минимизации определяется приближенным решением линейной системы с усеченной матрицей Гессе функционала (3), в которой отброшены члены, потенциально приводящие к потере положительной определенности. Итоговое смещение вдоль направления минимизации находится одномерным поиском.

2. ЗАДАНИЕ ФУНКЦИИ РАССТОЯНИЯ И ДРУГИХ ПАРАМЕТРОВ АДАПТАЦИИ

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

Методика построения восьмеричного дерева приведена в [15]. Начальная фоновая сетка разбивается до тех пор, пока минимальный размер ребра не будет меньше минимального сеточного размера желаемой адаптивной сетки.

Приближенные значения кривизны поверхности тела вычисляются на этапе подготовки фоновой сетки. Предполагается, что тело задано хорошо разрешенной триангулированной поверхностью. Для всех точек поверхностной сетки рассчитывается приближенное значение кривизны с использованием метода из работ [20], [21] (детальное описание дается в следующем разделе). Каждой вершине фонового восьмеричного дерева приписывается значение радиуса кривизны в ближайшей к данной вершине точке граничного контура/поверхности (см. фиг. 2).

3. МЕТОДИКА РАСЧЕТА ПОЛЕЙ КРИВИЗНЫ

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

4. ВЫБОР УПРАВЛЯЮЩИХ ПАРАМЕТРОВ

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

При задании сжатия в тангенциальном направлении требуется учет как кривизны границы тела, так и расстояния до медиальной оси. Использование только кривизны может приводить к резким переходам степени сжатия и, потенциально, к сеточным разрывам, особенно на внутренних углах (см. фиг. 4). Если расстояние до медиальной оси меньше, чем величина пограничного слоя в логическом пространстве, ${{\sigma }_{2}}$ на границе тела увеличивается, чтобы избежать “разрыва” сетки внутри тела. Увеличение ${{\sigma }_{2}}$ приводит к увеличению прообраза пограничного слоя в лагранжевых координатах, что приводит к увеличению числа вершин адаптивной сетки внутри области и около ее границы.

Рассмотрим алгоритм вычисления поля параметра ${{\sigma }_{1}}$, т.е. сжатия в нормальном направлении, где d – значение функции расстояния без знака до границы тела, h – длина ребра исходной сетки, оба параметра нормированы на L – величину по нормали адаптивной зоны:

Здесь $D$ – толщина переходного участка погранслойного слоя, $\delta $ – размер области по нормали, где сохраняется максимальное сгущение, $A$ – максимальная степень сжатия по нормали, ${{c}_{e}}$ – относительное приращение толщины двух соседних ячеек. Как правило, ${{c}_{e}}$ берется равным 0.2, чтобы не превышать рост сеточных элементов в нормальном направлении более, чем в 1.2 раза. Функция $d(\xi )$ на фиг. 5б показывает отображение из лагранжевых координат в эйлеровы в одномерном случае, а ${{\sigma }_{1}}(d) = ~\xi {\kern 1pt} '(d)$, что показано на фиг. 5а. Подробный вывод ${{\sigma }_{1}}(d)$ представлен в работе [10].

Для задания поля параметра ${{\sigma }_{2}}$ выполняются следующие вычисления.

Рассмотрим функцию для произвольного R:

где $L$ – размер по нормали адаптивной зоны, т.е. области, где узлы подвижны. Фиг. 6 иллюстрирует связь параметров из формулы (5), а именно: ${{\tilde {\sigma }}_{2}}(x,R) = {{l}_{\xi }}{\text{/}}{{l}_{x}}$.

Функция (5) имеет максимум не на границе. Чтобы это исправить, находим ${{d}_{{{\text{max}}}}}$ такое, что ${{\tilde {\sigma }}_{2}}({{d}_{{{\text{max}}}}}) = \max {{\tilde {\sigma }}_{2}}(d,~R)~$. Корректируем формулу (5):

Считаем ${{\sigma }_{{2c}}}~$и ${{\sigma }_{{2m}}}$: ${{\sigma }_{{2c}}} = {{\sigma }_{2}}(0,{{R}_{c}})$, где ${{R}_{c}}$ – радиус кривизны. Пусть ${{d}_{m}}$ есть расстояние от вершины контура до медиальной оси, а ${{\xi }_{m}} = \xi ({{d}_{m}})$ – расстояние от вершины до медиальной оси в лагранжевых переменных. Если ${{\xi }_{m}} < {{\xi }_{D}}$, то ${{\sigma }_{{2m}}} = \frac{{{{\xi }_{D}}}}{{{{\xi }_{m}}}}{{\sigma }_{{1,{\text{min}}}}}$. Если ${{\xi }_{m}} \geqslant {{\xi }_{D}}$, то ${{\sigma }_{{2m}}} = {{\sigma }_{{1,{\text{min}}}}}$. Из двух значений ${{\sigma }_{{2c}}}$ и ${{\sigma }_{{2m}}}~$выбираем максимальное. ${{\sigma }_{{{\text{2max}}}}} = \max ({{\sigma }_{{2c}}},{{\sigma }_{{2m}}})$ – значение на границе тела. Далее находим значение $R$ так, чтобы формула (6) была справедлива для найденного ${{\sigma }_{{2{\text{max}}}}}$. Иллюстрация вспомогательной функции ${{\tilde {\sigma }}_{2}}(d,R)$ и функции после корректировки ${{\sigma }_{2}}(d,R)$ для различных значений $R~$приведена на фиг. 7.

В трехмерной постановке добавляется третий параметр ${{\sigma }_{3}}$. Тогда построение полей ${{\sigma }_{2}}$ и ${{\sigma }_{3}}$ происходит аналогично, используя для задания $R$ значения главных кривизн.

6. ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ

ЗАКЛЮЧЕНИЕ

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

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

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