Метод Тейлора — Даффи для вычисления несобственных интегралов по треугольным панелям при решении граничных интегральных уравнений

Язык труда и переводы:
УДК:
519.64
Дата публикации:
28 ноября 2022, 21:25
Категория:
Математическое моделирование физических процессов и технических систем
Авторы
Иванова Юлия Витальевна
МГТУ им. Н.Э. Баумана
Хорошева Анна Александровна
МГТУ им. Н.Э. Баумана
Аннотация:
Рассмотрена реализация метода Тейлора — Даффи в контексте построения численной схемы повышенного порядка точности для решения граничного интегрального уравнения, возникающего в вихревых методах вычислительной гидродинамики при моделировании пространственного обтекания тел потоком несжимаемой среды. Использование алгоритма Тейлора — Даффи позволяет обеспечить разрывную кусочно-линейную аппроксимацию решения на кусочно-гладких обтекаемых поверхностях и реализацию схемы Галеркина для построения дискретного аналога граничного интегрального уравнения.
Ключевые слова:
вихревые методы, граничное интегральное уравнение, метод Тейлора — Даффи, схема Галеркина
Основной текст труда

Постановка задачи

Математические модели многих физических процессов сводятся к уравнению Лапласа или Гельмгольца. Одним из способов решения задач Дирихле и Неймана для этих уравнений является поиск решения в виде потенциала простого или двойного слоя [1]:


u({\vec {r}})=\int _{\partial \Omega }G({\vec {r}}-{\vec {\xi }})\mu ({\vec {\xi }})dS_{\xi },\quad u({\vec {r}})=\int _{\partial \Omega }{\frac {\partial G({\vec {r}}-{\vec {\xi }})}{\partial {\vec {n}}({\vec {\xi }})}}g({\vec {\xi }})dS_{\xi },\quad {\vec {r}}\in \Omega ,

где G — фундаментальное решение соответствующей задачи в двумерном или трехмерном случае. Интегральное представление решения через значения функции плотности потенциала на границе позволяет свести исходную задачу Дирихле к решению интегрального уравнения относительно плотности потенциала простого слоя: \int _{\partial \Omega }G({\vec {r}}-{\vec {\xi }})\mu ({\vec {\xi }})dS_{\xi }=f({\vec {r}}),

а задачу Неймана — к решению уравнения относительно плотности потенциала двойного слоя {\frac {\partial }{\partial {\vec {n}}({\vec {r}})}}\int _{\partial \Omega }{\frac {\partial G({\vec {r}}-{\vec {\xi }})}{\partial {\vec {n}}({\vec {\xi }})}}g({\vec {\xi }})dS_{\xi }=f({\vec {r}}).

Задачи вычислительной гидродинамики, связанные с моделированием обтекания профилей и тел, сводятся к задачам Неймана, а последнее уравнение может быть преобразовано к уравнению с гиперсингулярным (сильно-сингулярным) интегралом, понимаемым в смысле конечной части по Адамару: \int _{\partial \Omega }{\frac {\partial }{\partial {\vec {n}}({\vec {r}})}}{\frac {\partial G({\vec {r}}-{\vec {\xi }})}{\partial {\vec {n}}({\vec {\xi }})}}g({\vec {\xi }})dS_{\xi }=f({\vec {r}}).

В то же время вместо гиперсингулярного интегрального уравнения исходную задачу Неймана можно свести к уравнению относительно векторной величины градиента плотности потенциала двойного слоя [1, 2], которое имеет вид \int _{\partial \Omega }\nabla _{\xi }G({\vec {r}}-{\vec {\xi }})\times {\vec {\gamma }}({\vec {\xi }})dS_{\xi }-\alpha ({\vec {r}}){\vec {\gamma }}({\vec {r}})\times {\vec {n}}({\vec {r}})={\vec {f}}({\vec {r}}),\quad \quad \quad \quad \quad (1)

где {\vec {\gamma }}({\vec {r}})=-Grad(g({\vec {r}}))\times {\vec {n}}({\vec {r}}) имеет физический смысл интенсивности вихревого слоя на обтекаемой поверхности; оператор Grad(.)  означает поверхностный градиент, а интеграл является сингулярным (понимаемым в смысле главного значения по Коши) или несобственным.

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

Схема Галеркина

Приближенное решение уравнения (1) строится в виде линейной комбинации базисных функций, определенных на дискретизированной границе обтекаемой поверхности \partial \Omega , участки которой будем называть панелями и обозначать K_{i} . Дискретный аналог уравнения (1) получается из условия ортогонализации невязки (1) набору проекционных функций, которые в простейшем случае совпадают с базисными. 

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

Таким образом получаем, что коэффициенты линейной системы, выражающей дискретны аналог уравнения (1), даются двойными интегралами, вычисляемыми по паре панелей, вида \int _{K_{i}}\phi ^{(i)}({\vec {r}})dS_{r}\int _{K_{j}}G({\vec {r}}-{\vec {\xi }})\phi ^{(j)}({\vec {\xi }})dS_{\xi },\quad \quad \quad \quad \quad (2)

где \phi ^{(i)}({\vec {r}}) — функции указанного семейства.

При решении двумерных задач, когда граница \partial \Omega — замкнутая кривая, а панели K_{i} — прямолинейные отрезки, такие интегралы для кусочно-постоянных и кусочно-линейных базисных функций можно вычислить точно [3]. В трехмерных задачах, когда \partial \Omega — замкнутая поверхность, а панели K_{i} — плоские треугольники, — в [3] предложена методика вычисления интегралов для постоянных базисных функций; при этом внутренний интеграл в (2) можно вычислить аналитически, а внешний вычисляется численно. Если панели K_{i} и K_{j} не имеют общих точек, то использование квадратурных формул Гаусса позволяет легко решить эту задачу на практике, однако в противном случае внешний интеграл оказывается несобственным. В [3] разработана техника аддитивного выделения особенности, в соответствии с которой неограниченная составляющая интегрируется аналитически, а гладкая часть — численно, что в результате позволяет обеспечить достаточно высокую точность.

Повышение точности численного решения требует перехода от кусочно-постоянного к кусочно-линейному представлению решения, при этом обобщение методики [3] на случай линейных базисных функций представляется крайне трудоемким и едва ли рациональным. Вместо этого предлагается вычислять несобственные интегралы численно [4], производя так называемое преобразование Даффи, позволяющее в силу свойств ядра уравнения (1) сводить вычисление интегралов (2) к интегрированию ограниченных функций при достаточно произвольных базисных функциях \phi ^{(i)} .

Метод Тейлора — Даффи

Метод Тейлора — Даффи адаптирован для вычисления интегралов вида (2) по треугольным панелям K_{i} c линейными базисными функциями \phi ^{(i)} , считая, что каждой вершине панели соответствует одна базисная функция, обращающаяся в нуль в двух остальных вершинах. Такое введение базисных функций позволяет дать ясный геометрический смысл коэффициентам в разложении решения.

Применение метода Тейлора — Даффи необходимо в трех случаях:

  • панели K_{i} и K_{j} имеют общую вершину;
  • панели K_{i} и K_{j} имеют общую сторону;
  • панель K_{i}  совпадает с панелью  K_{j} .

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

Результаты

Для оценки точности метода Тейлора — Даффи и его вычислительной сложности было произведено сопоставление результатов расчетов для случая кусочно-постоянных базисных функций с результатами использования алгоритма [3], реализованного в коде Integrator [5] для сеток различного качества (рис. 1).

Рис. 1. Сетки высокого (слева), среднего (в центре) и низкого качества (справа)

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

Рис. 2. Среднее время интегрирования по паре панелей с общей вершиной (слева) и ребром (справа) методом Тейлора — Даффи (сплошные линии) и методом выделения особенности (штриховые линии)

На данный момент реализована схема решения гранитного интегрального уравнения пространственного обтекания тел на основе подхода Галеркина с применением метода Тейлора — Даффи для интегралов с кусочно-постоянными базисными функциями и схема метода Тейлора — Даффи для интегрирования по пространственному телу с использованием кусочно-линейных базисных функций.

Заключение

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

Литература
  1. Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент (в математической физике, аэродинамике, теории упругости и дифракции волн). Москва, Янус, 1995, 520 с.
  2. Kempka S.N., Glass M.W., Peery J.S., Strickland J.H., Ingber M.S. Accuracy considerations for implementing velocity boundary conditions in vorticity formulations. SANDIA report. SAND96-0583, UC-700. United States, 1996, 52 p. DOI: https://doi.org/10.2172/242701
  3. Марчевский И.К., Серафимова С.Р. Аналитическое и полуаналитическое вычисление интегралов от логарифмического и ньютоновского потенциала и их градиентов по прямолинейным отрезкам и треугольным панелям. Вычислительные методы и программирование, 2022, т. 23, № 2, с. 137–152.
  4. Тaylor D.J. Accurate and efficient numerical integration of weakly singular integrals in Galerkin EFIE solutions. IEEE Transactions on Antennas and Propagation, 2003, vol. 51, iss. 7, рp. 1630–1637. DOI: https://doi.org/10.1109/TAP.2003.813623
  5. Вычисление интегралов от фундаментального решения уравнения Лапласа и его градиента. Developer platform GitHub. URL: https://github.com/vortexmethods/integrator (дата обращения 09.09.2022).
Ваш браузер устарел и не обеспечивает полноценную и безопасную работу с сайтом.
Установите актуальную версию вашего браузера или одну из современных альтернатив.