Разработка и реализация Т-схем численного решения граничных интегральных уравнений в математических моделях вихревых методов вычислительной гидродинамики тема диссертации и автореферата по ВАК РФ 05.13.18, доктор наук Марчевский Илья Константинович

• Марчевский Илья Константинович
• доктор наук
• 2021

Оглавление диссертации доктор наук Марчевский Илья Константинович

Введение

Глава 1. Математические модели обтекания профиля, лежащие в основе вихревых методов

1.1. Математические модели обтекания профиля и тела, построенные на основе теории потенциала

1.2. Математическая модель обтекания профиля и тела, построенная на основе обобщенного разложения Гельмгольца

1.3. Математические модели пониженной размерности в плоских задачах обтекания профиля

1.4. Математические модели пониженной размерности в задачах пространственного обтекания тела

1.5. Точные решения модельных задач

1.5.1. Точные решения плоских модельных задач

1.5.2. Точные решения трехмерных модельных задач

1.5.3. Оценка точности приближенного решения

1.6. Выводы по главе

Глава 2. Расчетные схемы для решения плоских задач моделирования обтекания профилей, аппроксимируемых прямолинейными панелями

2.1. Классическая схема метода дискретных вихрей

2.2. Простейшая Т-схема для моделирования обтекания неподвижного профиля

2.3. Т-схема с кусочно-постоянным представлением решения

2.4. Построение дискретных аналогов интегральных уравнений

с помощью метода Галеркина

2.5. Т-схема с разрывным кусочно-линейным решением

2.6. Т-схемы с системами уравнений пониженной размерности

2.6.1. Экономичная схема с непрерывным решением

2.6.2. Экономичная схема с выделением разрыва решения

2.7. Т-схема с асимптотическим решением вблизи угловой точки

2.8. Выводы по главе

Глава 3. Расчетные схемы для решения плоских задач моделирования обтекания профилей, аппроксимируемых криволинейными панелями

3.1. Ограничения схем с прямолинейными панелями

3.2. Построение Т-схем с помощью метода Галеркина

3.3. Т-схема с кусочно-квадратичным представлением решения

3.4. Экономичные схемы

3.5. Коррекция решения для учета влияния вихрей

3.6. Выводы по главе

Глава 4. Расчетная схема для решения пространственных задач моделирования обтекания тел

4.1. Расчетная схема на основе Т-модели и подхода Галеркина

4.2. Дискретный аналог граничного интегрального уравнения

4.2.1. Вычисление интегралов в общем случае далеко расположенных панелей

4.2.2. Вычисление интегралов в случае общего ребра

4.2.3. Вычисление интегралов в случае общей вершины

4.3. Учет влияния завихренности в области течения

4.4. Восстановление плотности потенциала двойного слоя

4.5. Анализ точности разработанной Т-схемы

4.5.1. Анализ точности вычисления коэффициентов !ц

4.5.2. Анализ точности восстановления плотности потенциала двойного слоя

4.5.3. Анализ точности восстановления поля скоростей

4.6. Выводы по главе

Глава 5. Модификация алгоритма Барнса — Хата для быстрого приближенного матрично-векторного умножения

5.1. Проблемы непосредственного применения Т-схем

5.2. Алгоритм Барнса — Хата и его известная адаптация к решению двумерных задач вихревыми методами

5.3. Использование алгоритма Барнса — Хата для матрично-векторного умножения в двумерных задачах

5.4. Проблемы совместного использования алгоритма Барнса — Хата с Т-схемами повышенной точности

5.5. Модификация алгоритма Барнса — Хата на основе метода мультипольных разложений

5.6. Разложение скорости в двумерном случае

5.7. Разложение скорости в трехмерном случае

5.8. Выводы по главе

Глава 6. Алгоритмы итерационного решения линейных систем, возникающих в вихревых методах

6.1. Мультипольные моменты панелей

6.1.1. Мультипольные моменты прямолинейных панелей в плоском случае

6.1.2. Мультипольные моменты криволинейных панелей в плоском случае

6.1.3. Мультипольные моменты треугольных панелей в пространственном случае

6.2. Интегрирование локальных разложений по панелям

6.2.1. Вычисление интегралов вдоль прямолинейных панелей в плоском случае

6.2.2. Вычисление интегралов вдоль криволинейных панелей в плоском случае

6.2.3. Вычисление интегралов по треугольным панелям в пространственном случае

6.3. Итерационные методы решения линейных систем

6.4. Алгоритмы итерационного решения линейных систем в плоском случае

6.4.1. Предобуславливание системы для Т-схемы с кусочно-постоянным представлением решения

6.4.2. Умножение матрицы системы на вектор в Т-схеме с кусочно-постоянным представлением решения

6.4.3. Умножение матрицы системы на вектор в Т-схеме с кусочно-линейным представлением решения

6.4.4. Умножение матрицы системы на вектор в Т-схемах

при моделировании обтекания системы профилей

6.5. Алгоритмы итерационного решения линейных систем в пространственном случае

6.6. Выводы по главе

Глава 7. Программные реализации вихревых методов вычислительной гидродинамики

7.1. Программные комплексы VM2D и VM3D

7.1.1. Структура программных комплексов

7.1.2. Основные вычислительные блоки алгоритма

7.1.3. Технологии ускорения вычислений

7.2. Особенности реализации вихревых методов при использовании Т-схем

7.2.1. Преобразование завихренности в форме распределенного вихревого слоя в вихревые элементы

7.2.2. Определение скорости среды вблизи обтекаемой поверхности

7.3. Расчет присоединенных масс профилей и тел

7.3.1. Методика расчета присоединенных масс

7.3.2. Расчет присоединенных масс профиля

7.3.3. Расчет присоединенных масс трехосного эллипсоида

7.4. Примеры моделирования обтекания профилей и расчет действующих на них нагрузок

7.4.1. Моделирование течения Блазиуса

7.4.2. Моделирование течения в канале с обратным уступом

7.4.3. Определение положения точки отрыва потока при обтекании кругового цилиндра

7.4.4. Моделирование обтекания крыловых профилей

7.4.5. Моделирование нестационарного обтекания движущегося цилиндра

7.4.6. Моделирование обтекания системы близко расположенных профилей

7.5. Примеры моделирования обтекания тел и расчет действующих на них нагрузок

7.5.1. Моделирование обтекания сферы

7.5.2. Моделирование обтекания крыла конечного размаха

7.6. Выводы по главе

Общие выводы по работе и заключение

Список литературы

452

Введение диссертации (часть автореферата) на тему «Разработка и реализация Т-схем численного решения граничных интегральных уравнений в математических моделях вихревых методов вычислительной гидродинамики»

Введение

Актуальность темы исследования. Вихревые методы вычислительной гидродинамики являются, по-видимому, одними из первых численных методов, которые стали применяться при решении актуальных практических задач. Первые их модификации и реализации появились еще в 50-е годы XX века, практически одновременно с появлением ЭВМ. Они явились, по существу, следующим шагом в эволюции полуаналитических подходов, развивавшихся пионерами аэродинамики — О. Лилиен-талем, Н. Е. Жуковским, С. А. Чаплыгиным, Л. Прандтлем, М. В. Куттой, Т. Карманом, Л. Розенхедом и многими другими отечественными и зарубежными учеными, результаты фундаментальных исследований которых дали возможность превратить конструирование летательных аппаратов и выполнение полетов из наивных опытов по подражанию полету птиц в строгие научные дисциплины, а интуицию первых конструкторов и авиаторов заменить научно-обоснованными принципами конструирования и аэродинамического расчета.

Интересно проследить «научную биографию» основоположников отечественной научной школы вихревых методов, которая восходит напрямую к Николаю Егоровичу Жуковскому (1847-1921 гг.).

Сам Н. Е. Жуковский, несомненно, принадлежит к числу блестящих русских ученых, работы которых являются поистине пионерскими и во многом определили образ современного мира. Ему покорились те задачи, над которыми работали и которые не смогли решить Ньютон и Леонардо да Винчи: он смог построить теорию подъемной силы крыла и тянущей силы воздушного винта. Именно Н. Е. Жуковский в 1898 году, за 5 лет до полета «Флайера» братьев Райт, сказал:

Человек не имеет крыльев и по отношению веса своего тела к весу мускулов в 72 раза слабее птицы... Но я думаю, что он полетит, опираясь не на силу своих мускулов, а на силу своего разума.

Прошло совсем немного времени, и в 1906 году он сделал решающий шаг, чтобы самому осуществить мечту человечества о полетах: выпустил работу «О присоединенных вихрях» [50], которая, несмотря на скромное название, открыла, по существу, новую эпоху: именно в ней была блестяще изложена новая теория, позволившая впервые объяснить возникновение подъемной силы крыла и перейти к конструированию летательных аппаратов на основе строгого расчета. Фактически, им установлена возможность замены обтекаемого профиля системой вихрей, связанных с профилем и поэтому названных «присоединенными», которые оказывают на течение среды такое же влияние, как и обтекаемый профиль. Ему принадлежит формула для расчета подъемной силы профиля крыла (известная как теорема Н. Е. Жуковского), которая приведена в упомянутой работе и устанавливает именно вихревую природу подъемной силы.

Сложнейшие проблемы механики он, магистр чистой математики, профессор Московского университета и Императорского Московского технического училища (ныне — Московского государственного технического университета имени Н.Э. Баумана), смог решить благодаря изящному применению современных ему математических методов, и этим определил для современников и последователей своеобразный идеал ученого-механика и одновременно математика-прикладника, актуальный, несомненно, и сегодня. Едва ли он мог представлять себе, какие возможности откроют исследователям ЭВМ через полвека, через век, и какое место в науке и технике будут занимать методы расчета и математического моделирования, но тот факт, что многие современные методы расчетов в аэродинамике напрямую восходят к работам Н. Е. Жуковского, лишний раз подтверждает масштаб его научного гения.

Благодаря инициативе Н. Е. Жуковского в нашей стране были созданы такие центры авиационной науки, как Центральный аэрогидродинамический институт (ЦАГИ) и Военно-воздушная инженерная акаде-

мия (ВВИА), получившие всемирную известность и ныне носящие его имя. Он создал замечательную школу, воспитал целую плеяду учеников, ставших блестящими учеными, инженерами, авиаконструкторами.

Прежде всего, следует упомянуть работы Сергея Алексеевича Чаплыгина (1869-1942 гг.), прямого ученика Н.Е. Жуковского, чье имя увековечено, в частности, в сформулированном им совместно со своим учителем постулате (условие Чаплыгина — Жуковского, или Чаплыгина — Жуковского — Кутты, поскольку к аналогичному результату пришел также и немецкий математик М. В. Кутта) [123], который позволяет определять циркуляцию поля скоростей, а следовательно — и подъемную силу при обтекании профиля с острой кромкой. Данное условие было затем положено в основу многих расчетных схем вихревых методов.

Одним из учеников С. А. Чаплыгина был Николай Евграфович Кочин (1901-1944 гг.), знаменитый своими работами как в области гидродинамики и газовой динамики, так и развитием математического аппарата векторного и тензорного исчисления, который, восходя к классическим результатам Г. Гельмгольца, А. Фридмана и других исследователей, лежит в основе многих современных математических моделей.

Владимир Васильевич Голубев (1884-1954 гг.), хотя и не являлся непосредственным учеником Н. Е. Жуковского и С. А. Чаплыгина, но, получив блестящее математическое образование под руководством Д. Ф. Егорова и Н. Н. Лузина, уже с середины 1920-х годов активно занялся проблемами аэромеханики после знакомства с работами Жуковского и Чаплыгина. Ему, в частности, принадлежат работы по теории крыла с механизацией, теории машущего полета, а также фундаментальный труд «Лекции по теории крыла». С 1932 до 1954 года В. В. Голубев возглавлял кафедру высшей математики Военно-воздушной инженерной академии им. проф. Н. Е. Жуковского, являясь одновременно создателем и первым заведующим кафедрой аэромеханики, а также первым деканом механико-математического факультета Московского университета.

Учеником В. В. Голубева являлся Сергей Михайлович Белоцерков-ский (1920-2000 гг.), закончивший академию в 1945 году и занявшийся в ней активной научно-исследовательской деятельностью. Впоследствии он возглавил кафедру аэродинамики и был заместителем начальника академии. Именно С.М. Белоцерковского с полным основанием можно считать основоположником отечественной научной школы вихревых методов вычислительной аэрогидродинамики. Им в 1950-х годах для решения задач аэродинамики был разработан метод дискретных вихрей, в основе которого лежит сведение задачи расчета обтекания контура или поверхности к решению сингулярного граничного интегрального уравнения [11]. Теория сингулярных интегральных уравнений восходит к трудам Д. Гильберта и А. Пуанкаре, а одной из первых работ, в которой задача обтекания тонкого профиля была сведена к сингулярному уравнению, была работа М. А. Лаврентьева [70], в ней же содержались некоторые идеи относительно методов приближенного численного решения таких уравнений.

Метод дискретных вихрей, предложенный С. М. Белоцерковским, можно рассматривать как реализацию метода коллокаций при численном решении сингулярных интегральных уравнений совместно с применением квадратурных формул, аналогичных, по существу, методу центральных прямоугольников. Не ставя перед собой задачу обзора колоссального объема исследований, проведенных С. М. Белоцерковским и при его участии, упомянем среди множества лишь некоторых его ближайших учеников и соавторов — М. И. Ништа, Б. Г. Скрипача, В. Г. Табачникова, В.А. Апаринова, В.Н. Котовского, Р. М. Федорова, А. И. Желанникова, Б. С. Крицкого, В. В. Гуляева. Ими развиты многочисленные модификации метода дискретных вихрей и замкнутых вихревых рамок — бессеточных чисто лагранжевых методов моделирования плоских и пространственных течений [12-14,91,96,102]. Внедрение этих методов в практику расчетов позволило решить актуальные на то время задачи расчета характеристик новых типов летательных аппаратов и новых конструкций, к примеру, решетчатого крыла. Лагранжев характер вихревых ме-

тодов определяется фундаментальными свойствами поля завихренности, которые были установлены Гельмгольцем [93, 169] и означают «вморо-женность» вихревых линий в идеальную несжимаемую среду.

Применение метода дискретных вихрей и выполнение расчетов на ЭВМ первых поколений в 1950-60-е годы стало возможным благодаря исключительно малой, по современным меркам, трудоемкости вычислительного алгоритма. В его основу положены наиболее простая из возможных моделей среды — идеальная несжимаемая среда и некоторые гипотезы, в том числе аналогичные постулату Чаплыгина — Жуковского. Указанные ограничения сужают область применимости метода, однако решение актуальных на то время задач стало возможным с приемлемой для практических целей точностью.

Вопросами математического обоснования вихревых методов занялся в Академии им. Н. Е. Жуковского на кафедре Высшей математики Иван Кузьмич Лифанов (1942-2016 гг.), который впоследствии возглавил эту кафедру и стал крупным специалистом в области сингулярных интегральных уравнений. Ему принадлежат фундаментальные результаты, позволившие вскрыть основы и дать строгое математическое описание моделям и методам, лежащим в основе вихревых методов [12,21,76]. Развитый им математический аппарат и методы численного решения сингулярных интегральных уравнений, являясь развитием работ Н.И. Му-схелишвили [95] и Ф. Д. Гахова [27], с успехом применяются сегодня при решении широкого диапазона задач математической физики: от дифракции и рассеяния акустических и электромагнитных волн до фильтрации жидкости в пористых средах и теории упругости. Вместе с Юрием Владимировичем Ганделем (1934-2017 гг.) им в 1983 году был организован всесоюзный, а впоследствии — международный симпозиум «Метод дискретных особенностей в задачах математической физики», ставший наряду с научным семинаром С. М. Белоцерковского (позднее — имени С.М. Белоцерковского), действующим с 1959 года, фактически центральной площадкой для обсуждения как математических аспектов, так и приложений соответствующих методов. Заслуживают внимания в этом

направлении многочисленные работы Л.Н. Полтавского, Г. М. Вайник-ко, А. В. Сетухи, В. Ф. Пивня, Д. Г. Саникидзе, Ш. С. Хубежты и других исследователей. Из зарубежных работ в области обоснования вихревых методов следует отметить работу [140] и другие работы этих же авторов.

Не снижая значимости работ зарубежных исследователей, восходящих в значительной степени к пионерской работе автора R. Rosenhead [258], обзоры которых по состоянию на 1980-е годы содержатся в работах авторов A. Leonard [212] и T. Sarpkaya [264], а более свежих работ — в статье [137], монографиях авторов R. Lewis [216], G.-H. Cottet и P. Koumoutsakos [147], E. Branlard [141] и обзоре авторов C. Mimeau и I. Mortazavi [272], отметим, что они в основном посвящены моделированию эволюции завихренности в области течения. Следует отметить, что данные вопросы изучены к настоящему времени достаточно подробно, накоплен также и достаточный опыт применения различных численных схем и вычислительных алгоритмов.

Так, в начале 1970-х годов А. Чориным [145] был предложен так называемый метод случайных блужданий (random walk), основанный на возможности описания процесса диффузии завихренности в вязкой жидкости с помощью вероятностных моделей, восходящих к уравнению Фок-кера — Планка — Колмогорова, которому также подчиняются характеристики процесса броуновского движения частиц в соответствии с теорией А. Эйнштейна [157]. Применение такого подхода позволило использовать вихревые методы не только в рамках простейшей модели среды — идеальной (невязкой) несжимаемой жидкости, но и учитывать влияние вязкости среды. При этом интенсивности вихревых частиц в данном методе остаются неизменными, а к конвективной скорости их движения (т. е. скорости среды) прибавляется случайная величина с нулевым математическим ожиданием и дисперсией, пропорциональной вязкости среды. При достаточно большом количестве вихревых частиц и малом шаге расчета по времени данный метод позволяет обеспечивать высокую точность моделирования течения вязкой среды. Пример использования данного подхода в рамках современных реализаций вихревых методов представлен в работе [162].

Позднее P. Degond и S. Mas-Gallic был разработан метод обмена интенсивностями (particle strength exchange, PSE) [148], в котором интенсивности вихревых частиц не остаются постоянными, а перераспределяются на каждом шаге расчета. Оба указанных метода получили широкое распространение [147,244,257,267,270,290] и активно применяются до настоящего времени.

Определенную популярность получил так называемый метод расширения ядра (core spreading), описанный в выше упоминавшемся обзоре [212], суть которого сводится к представлению каждой вихревой частицы гауссовым распределением завихренности с изменяющейся со временем дисперсией, однако его первоначальная реализация оказалась не вполне корректной, на что указал C. Greengard [163]. Корректная модификация данного метода была разработана L. Rossi [259]; ядро вихря расширяется лишь до определенного предела, после чего вихрь разделяется на несколько вихрей меньшего размера.

Большой интерес представляет метод диффузионной скорости (diffusion velocity), предложенный Y. Ogami и T. Akamatsu [248], который, будучи в некоторой степени аналогичным методу случайных блужданий, является при этом полностью детерминированным. Согласно нему, вихревые частицы движутся в плоских течениях вязкой несжимаемой жидкости вдоль поля скоростей, представляющего собой суперпозицию конвективной скорости (скорости среды) и диффузионной скорости, пропорциональной вязкости среды и характеризующей неравномерность распределения завихренности в области течения. Отметим, что исследования в этом направлении восходят к фундаментальным работам А. Фридмана [93,127], в связи с чем диффузионную скорость в [120] предложено называть «скоростью Фридмана».

Метод диффузионной скорости получил существенное развитие в работах Г. Я. Дынниковой [3,42,44], в которых предложена оригинальная

методика вычисления диффузионных скоростей вихревых частиц, позволяющая корректно учитывать различную их концентрацию в различных частях области течения, которая может также существенно изменяться во времени. Данный метод положен в основу метода вязких вихревых доменов (ВВД) [3,35,41,46] — численного метода моделирования плоских течений вязкой несжимаемой среды, который оказался весьма эффективным инструментом решения широкого класса задач, в том числе в сопряженной постановке. Следует также отметить фундаментальные работы Г. Я. Дынниковой, в которых получен аналог классических интегралов Бернулли и Коши — Лагранжа [40], предложены методы расчета гидродинамических нагрузок, действующих на обтекаемые тела [154,155], а также установлена структура гидродинамических сил, действующих на тело, движущееся в жидкости [45]. Указанные работы имеют колоссальное значение для развития вихревых методов применительно к моделированию как плоских, так и пространственных течений несжимаемой среды. Метод диффузионной скорости получил дальнейшее развитие в работах А. В. Сетухи [118], где он фактически получил строгое математическое обоснование и был перенесен (при некоторых упрощающих предположениях) на случай моделирования пространственных течений.

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

Заслуживает внимания метод изолированных точечных вихревых частиц — вортонов, предложенный в работе Е. А. Новикова [97], однако главный его недостаток заключается в несоленоидальности поля завихренности, представленного изолированными вихревыми частицами, если только данные вихревые элементы не располагаются в пространстве строго определенным образом, образуя замкнутые вихревые структуры; однако даже если это и так, то в процессе движения такое расположение будет нарушено вместе с условием бездивергентности завихренности. Аналогичные модели обсуждались также в работе [130].

Обобщением модели точечного вортона можно считать модель вихревого отрезка. Вихревые отрезки рассматривались еще в рамках метода вихревых рамок в качестве составляющих самих рамок, однако теперь они могут выступать в качестве изолированных вихревых элементов. Данный подход использовался в работах А. В. Сетухи, Г. А. Щеглова, И. К. Марчевского [6,16,30,89,230], однако вычислительные эксперименты показали, что удовлетворительных результатов можно добиться лишь применяя специальные процедуры коррекции, позволяющие «уменьшить» степень несоленоидальности поля завихренности в расчете [231].

Решить проблему несоленоидальности поля завихренности удается путем представления завихренности в области течения в виде замкнутых вихревых структур — вихревых петель. Такой подход описан в работах [283,289], при этом следует отметить, что первоначально он применялся для создания анимаций и визуализаций течений жидкостей и газов. Тем не менее, его удалось успешно применить и для решения задач пространственного моделирования обтекания тел и расчета действующих на них гидродинамических нагрузок; подобные модификации разработаны С. А. Дергачевым, И. К. Марчевским и Г. А. Щегловым [37,38,149].

Представляет интерес предложенный Г. Я. Дынниковой метод вихревых диполей [41,47], который хотя и не относится непосредственно к вихревым методам (первичной расчетной величиной в нем является

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

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

В случае идеальной жидкости обеспечение непротекания достигается расположением на обтекаемой поверхности вихревого слоя (или двойного слоя; данные понятия тесно связаны друг с другом, подробное обсуждение данной взаимосвязи имеется также в [76]). Задача определения интенсивности вихревого (или двойного) слоя, математически аналогичная внешней или внутренней задаче Неймана (для неограниченных и ограниченных течений соответственно), традиционно решается методами теории потенциала [21,69,76,122,173] и сводится к решению сингулярного или гиперсингулярного граничного интегрального уравнения. При этом течение среды остается безвихревым, за исключением области вихревого следа, формирующегося в виде тонких вихревых пелен, сходящих с острых кромок (в соответствии с упоминавшимся выше условием Чаплыгина — Жуковского — Кутты) или каких-либо еще априори заданных точек.

Собственно, И. К. Лифановым предложены и обоснованы численные методы решения таких уравнений, получившие название квадратур-

ных формул типа дискретных вихрей [12,21,76]. Применительно к решению пространственных задач методом вихревых рамок обоснованием метода также занимается А. В. Сетуха с коллегами [57,105,106,119]; большое количество расчетных схем повышенной точности в рамках данного подхода построено его учениками [36,99,110].

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

При моделировании течений вязкой жидкости на обтекаемой поверхности выполняется условие прилипания. Это вызывает генерацию завихренности вблизи обтекаемой поверхности и ее поток в область течения; понятие потока завихренности введено в [218], связанные с этим процессы обсуждаются в монографии [20]. Завихренность в этом случае нельзя считать присоединенной, поскольку она фактически уже является частью области течения. Течение, соответственно, перестает быть безвихревым (и тем более — потенциальным); вблизи обтекаемой поверхности образуется некоторое распределение завихренности. В этой связи следует отметить, что в численных методах — не важно, сеточных или бессеточных, — в которых первичной расчетной величиной является завихренность, сформулировать граничное условие на обтекаемой поверхности в терминах завихренности не удается, а имеющиеся подходы носят приближенный характер. Поэтому граничное условие ставится для поля скоростей, а его удовлетворение обеспечивается потоком завихренности, интенсивность которого определяется из решения некоторого граничного интегрального уравнения. Подобные уравнения записаны в уже упоминавшихся работах Г. Я. Дынниковой [41,45] и А. В. Сетухи [118]. С учетом оценок, сделанных в [20], можно перейти к задаче определения интенсивности тонкого вихревого слоя, который образуется за малый

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

В то же время следует отметить, что использование разработанных И. К. Лифановым численных схем решения сингулярных интегральных уравнений совместно с современными модификациями вихревых методов может приводить к получению результатов с намного меньшей точностью, чем в методе дискретных вихрей (хотя даже там, как установлено в [76], они обеспечивают порядок точности не выше первого). Это связано с тем обстоятельством, что в области течения вблизи обтекаемой поверхности теперь имеется завихренность, моделируемая большим количеством вихревых элементов. Их влияние на решение интегрального уравнения выражается в структуре функции правой части, которая теперь может иметь значительный поверхностный градиент. Кроме того, требуется повышенное качество восстановления поля скоростей в пристеночной области, поскольку от этого напрямую зависит правильность моделирования эволюции находящейся там завихренности, а она, в свою очередь, влияет на величину потока завихренности и, таким образом, на интенсивность ее образования вблизи обтекаемой поверхности. О систематических исследованиях в этом направлении не известно, однако факт невысокой точности существующего подхода описан в [174] и подтвержден вычислительными экспериментами в модельных задачах в диссертации [94]. В той же работе [94] установлено, что именно невысокая точность решения граничного интегрального уравнения и, как следствие,

погрешность определения интенсивности вихревого слоя является наиболее существенным источником погрешности в методе вязких вихревых доменов, тогда как повышение точности численного интегрирования уравнений движения вихревых элементов [84, 94] оказывает существенно меньшее влияние на результат моделирования обтекания профилей. Тот же вывод можно сделать из анализа представленных в работе [37] результатов решения трехмерных задач с помощью метода вихревых петель.

Альтернативой, предложенной в [174], является рассмотрение граничного интегрального уравнения 2-го рода относительно интенсивности вихревого слоя. Свойства этого уравнения в указанной работе не обсуждаются, однако можно показать, что его ядро является ограниченным или, по крайней мере, абсолютно интегрируемым. Данную математическую модель назовем T-моделью. Отметим, что идентичную структуру имеет интегральное уравнение относительно потока завихренности, полученное в [118]. На возможность рассмотрения такой математической модели указано в обзоре [264], а также в монографиях [76,147], однако ни о каких дальнейших исследованиях в этом направлении не известно, если не считать работы [158], не получившей, впрочем, систематического развития, в которой фактически был использован указанный подход.

Список литературы диссертационного исследования доктор наук Марчевский Илья Константинович, 2021 год

r - —

Поле скоростей, индуцируемое свободным вихревым слоем в области течения

Vy(р) = jf (y(£) х Q(P - —)) dSe, r e F,

должно быть потенциальным, поэтому должно выполняться условие [76]

Div Y(r) = 0, r e K,

уже упоминавшееся в разделе 1.1, где оператор Div означает поверхностную (двумерную) дивергенцию соответствующего распределения. Последнее означает, что интенсивность вихревого слоя Y(r) может быть выражена через поверхностный градиент потенциала двойного слоя

Y(r) = — Gradg(r) х n(r), r e K,

относительно которого можно записать эквивалентное гиперсингулярное интегральное уравнение (1.14), которое в трехмерных задачах при G(r, —) = ^—— принимает вид

/ (n(r) • n(—)) |r - —|2 - 3(n(r) • (r - —))(n(—) • (r - —))

fK-4^-g(—) dS = fn(r1

r e K, (1.29)

и интеграл в котором следует понимать в смысле конечной части по Адамару [2,76].

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

Vт(г) = V(г) - (V(г) ■ п(г))п(г) = п(г) X (V(г) х п(г)).

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

Таким образом, интегральное уравнение (1.21) в проекции на касательную плоскость можно записать в виде

Цп(г) X ( 7(£|Х -г ^0 X п(г))) dSe - а (г )(7(г) х п(г)) = / т (г),

г е К, (1.30)

или

£ Р т (г, 0) dSe - а (г )(7(г) X п(г)) = / т (г), г е К,

'к

где правая часть имеет вид

/ т (г) = п(г) X ((а(г )и К (г) - V (П(£) х Я(г - 0) dVí -

-£ (т*^) X Я (г - £)) dSe - ^ ц*{{(£Жг - £)) dSe) X п(г )).

Характерной чертой данного интегрального уравнения является интегрируемость его ядра: хотя оно и является неограниченным при |г - ^ 0, поверхностный интеграл в нем сходится как несобственный, поскольку порядок особенности в подынтегральном выражении снижается на единицу в результате векторного умножения на орт нормали.

Рассмотрим для примера случай, когда в окрестности некоторой точки г о е К обтекаемая поверхность является гладкой и принадлежит классу С3.

Рассмотрим локальную декартову систему координат Охуг, оси Ох и Оу которой лежат в касательной плоскости, а ось Ог направлена по нормали к поверхности (Рис. 1.4).

Рис. 1.4. Участок гладкой обтекаемой поверхности, касательная плоскость и локальная система координат

Тогда можно считать, что форма поверхности в некоторой окрестности точки г о может быть задача уравнением г = / (х, у), при этом

/(о, о) = о, /х(о, о) = /у(о, о) = о

и, следовательно, с точностью до малых третьего порядка отклонение поверхности от касательной плоскости в окрестности точки касания определяется второй квадратичной формой поверхности [103], матрицу которой обозначим за 12, а ее компоненты — за Ь, М, N:

г = 2(Ьх2 + 2Мху + Ny2) = 1 а%а,

где а = п(го) х ((го — 0 х п(го)) — компонента вектора (го — 0, лежащая в касательной плоскости.

Поскольку вектор интенсивности вихревого слоя 7(0 лежит в касательной плоскости, его третья компонента есть величина порядка |а(0|2, тогда в окрестности точки го справедливо приближенное равенство

а^ а

р т ^ 003Пг о) х 7(г о)),

которое можно проинтегрировать по области, проецирующейся на круг радиусом е на касательной плоскости; в результате получаем

Се £

JK P т dSe « 8 (L + N К 0) х Y(r 0)) =8 tr I2 (n(r о) x 7(r 0)),

где tr I2 обозначает первый инвариант (след) матрицы.

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

Div Y(r) = 0, r Е K, из которого, в частности, следует, что

jf7(0 dSe = 0.

В этой связи использование при решении пространственных задач N-модели и сведение получающегося уравнения (1.11) к гиперсингулярному скалярному интегральному уравнению первого рода относительно плотности потенциала двойного слоя (1.14), которое при d = 3 принимает вид (1.29), представляется более логичным. Если предположить, что уравнение (1.29) тем или иным способом решено и поверхностное распределение плотности потенциала двойного слоя g(r), r Е K, известно, то интенсивность вихревого слоя

7(r) = — Gradg(r) x n(r), r Е K,

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

Недостатком такого подхода является сравнительно низкая, как отмечено в [21], точность решения гиперсингулярного граничного интегрального уравнения, соответствующего N-модели, особенно для обтекаемых тел сложной формы, поскольку квадратурные формулы, позволяющие выделять конечные части гиперсингулярных интегралов, изначально имеющие невысокий порядок точности [21,76], оказываются весьма «чувствительными» к качеству поверхностной сетки. Ранее во Введении было упомянуто об исследованиях, в результате которых были построены алгоритмы повышенной точности, позволившие в некоторой мере снизить остроту проблемы неравномерных неструктурированных поверхностных сеток невысокого качества, однако решить ее принципиально подобным образом, по-видимому, не удастся.

В то же время представляется, что использование Т-модели (1.30) позволит добиться значительного повышения точности определения интенсивности вихревого слоя 7(г), также как и плотности потенциала двойного слоя д(г), г е К. Именно факт абсолютной интегрируемости ядра открывает перспективу разработки вычислительных алгоритмов, которые способны обеспечить высокую точность решения даже при сравнительно грубой аппроксимации обтекаемой поверхности и низком качестве поверхностных сеток. При этом следует помнить, что в силу неизбежной погрешности аппроксимации условие бездивергентности численного решения (1.10) в практических расчетах будет нарушено.

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

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

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

• для N-модели — гиперсингулярного скалярного граничного интегрального уравнения относительно плотности потенциала двойного слоя с последующим, при необходимости, вычислением интенсивности вихревого слоя;

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

1.5. Точные решения модельных задач

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

1.5.1. Точные решения плоских модельных задач

Наиболее простым способом получения точных решений плоских модельных задач является метод конформных отображений [71,75,93]. Его можно применить для профилей простых форм: эллипсов, крыловых профилей Жуковского, а также некоторых других.

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

Для построения точного решения для интенсивности 7(г) вихревого слоя на границе профиля (а если требуется — то и для скорости среды V(р) в области течения) сначала запишем выражение для комплексного потенциала течения при моделировании обтекания кругового профиля радиусом Я в присутствии точечных вихрей в потоке.

Будем считать, что центр окружности С радиусом Я располагается в начале координат комплексной плоскости (£), Wc» — комплексная скорость набегающего потока, т. е. ее действительная и мнимая часть равны соответственно горизонтальной и вертикальной компонентам истинной скорости набегающего потока. Примем, что точечные вихри в количестве N0 штук, имеющие циркуляции Гт, расположены в точках = рте1°т комплексной плоскости, т = 1, ..., N1.

Чтобы учесть вклад влияния т-го вихря в комплексный потенциал течения, следует построить систему отраженных вихрей (Рис. 1.5), ко-

торая состоит из вихря с циркуляцией противоположного знака (—Гт), располагаемого в симметричной по отношению к окружности точке

г?2

4« = ■ е">■ = -е">■, (1.31)

рт

и вихря с циркуляцией Гт, помещаемого в начало координат — в центр окружности.

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

Рис. 1.5. Круговой профиль, точечный вихрь в области течения и система отраженных вихрей

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

выполнение перечисленных выше свойств:

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

Для профиля произвольной формы, ограниченного контуром К, комплексный потенциал течения можно построить, если известны конформные отображения 2 = £(£) и £ = £(2), переводящие внешность окружности С на внешность контура К и наоборот соответственно (Рис. 1.6), удовлетворяющие условию на бесконечности

Рис. 1.6. Круговой профиль С на плоскости (£) и профиль произвольной формы, ограниченный контуром К, на плоскости (2) и конформные отображения внешностей этих профилей

Чтобы записать выражение для комплексного потенциала течения вне кругового профиля С на плоскости (£), который будет соответствовать комплексному потенциалу обтекания исходного контура К на плоскости (2), необходимо отобразить точечные вихри из плоскости (2) на

('(£) ^ сопб1 при £ ^ то.

О

г

плоскость %

плоскость г

плоскость (£), и затем построить систему отраженных вихрей. В результате для точек 2, лежащих во внешности контура К, можно записать

Р (О = Р (&)) = / (2) = WU(z) + + 2П1п 1Г +

+ £ (,п^ ^2) - ,пС(2) - ^^ + |пМ) . (,.32)

■=1 п V /

Здесь 2т, ■ = 1, ..., N — положения точечных вихрей на комплексной плоскости 2 в области течения; комплексная скорость Ж» находится как

Ж» = К»- С(»),

где V» — комплексная скорость набегающего потока в исходной задаче (при обтекании контура К).

Для эллиптического профиля и крылового профиля Жуковского необходимые конформные отображения известны и строятся на основе функции Жуковского [71,75,93]:

С(а = ! (с + н + О+Н) ' ^) = 2 - Н + - а2, (1.33)

где параметры а и Н, а также радиус Я окружности С могут быть определены по геометрическим параметрам исходного профиля:

• для эллиптического профиля

а

= у/а2 - Ь2, Я = а1+ Ь1, Н = 0,

где а1 и Ь1 — полуоси исходного эллипса; для профиля Жуковского

Я = й + уОТ^, н = -ай + Ш, (1.34)

л/а2 + Н2

где параметры а, й и Н связаны с длиной, толщиной и кривизной крылового профиля соответственно.

Заметим, что для введенного таким образом отображения

С(£) ^ 2 при £ ^ то,

поэтому Што = 2Уто.

Далее на окружности С в плоскости (£) удобно ввести параметризацию с помощью параметра г, тогда

£ = г е [0, 2п),

и для комплексного потенциала на границе профиля К можно записать

ф(г) = ^(Я е^-^) = ШтоЯ е1(г-{р) + ШтоЯ е-1(г-{р) +

/ _ ш / ч Л/оГ я ;(—) _ еюш

+ Г + У гЛ + у Гш 1п —-, г е [0, 2п),

где рше1°ш = £(2ш) — положение ш-го точечного вихря после его отображения на плоскость (£); рЩ?) = Я2/рш.

Величину ^ удобно выбрать таким образом, чтобы значение параметра г = 0 соответствовало точке на большой полуоси эллипса (тогда ^ = 0) или острой кромке на профиле Жуковского (в этом случае следует принять ^ = аг^ ^).

На исходном контуре К значению параметра г, которое отвечает некоторой точке окружности С на плоскости (£), соответствует точка

1 / _ а2

2

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

2(0 = с (е^) = 2( х(0 + ^ )' х(г) = Я е^ + н. (1.35)

комплексное число единичного модуля, задающего направление касательной:

7 (г) = V(0 • т(г) = —— • —— = ——.

г'(г) \г'(0\ \г'(0\

Окончательно можно получить формулу для точного решения уравнения, выражающего интенсивность вихревого слоя на контуре [201]:

*(А 2\у»\ 8т(^ + в - г) + (г + г(р)/М) (1 т (г) =-^---, (1.36)

а

1 а

(Яе1(—)+ И)'

где величина ^(г) определяет влияние точечных вихрей в области течения:

г (г) = £ н С08(г - - у,1)~ ^ г..

да=1 ^ соБ(г - от - - - + р^)

Величины pw и Qw вычисляются по формулам

Pw = |C(2w )|, #w = arg Z (2w),

при этом следует выбрать ту ветвь функции квадратного корня из комплексного числа в функции Z(^), которая обеспечивает выполнение условия pw > R. Величины pR находятся по правилу (1.31).

Отметим, что знаменатель в выражении (1.36) можно выразить явно через параметры профиля:

для эллиптического профиля

2

О

1 О

(Rei(—)+ H )2

2 у/о2 sin21 + b2 cos21

oi + bi

для профиля Жуковского

2 О2

(Rei(i-^)+ H)'

2

2 2 (d sin ф + h + R sin(t — ф)) + (R cos(t — ф — d cos ф

x

| ((d sin ф+R sin(t — ^))(d sin ф+h+R sin(t — ф)) + dR(cos t — 1)) +

/ \211/2 + (R(sin(t — ф) + sin ф) (R cos(t — ф) — d cos ф) J > .

Величина циркуляции поля скоростей Г вокруг профиля, равная суммарной завихренности, содержащейся в вихревом слое, для эллиптического профиля может быть выбрана произвольно; в рассматриваемых тестовых задачах полагается Г = 0 независимо от угла атаки профиля.

Для крылового профиля Жуковского величину циркуляции Г можно найти из условия Чаплыгина — Жуковского об ограниченности скорости на острой кромке:

Na

Г = -2nR| sin(a + - У

R cos(#w + <р) - pW

(r)

(r)

W=1 R COS(^w + <p) - 1 (pw + pw где а = arg — угол атаки профиля.

)

-г

w

(1.37)

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

На Рис. 1.7, а, изображены линии тока при обтекании эллиптического профиля с полуосями а1 = 4 и Ь1 = 1 набегающим потоком, имеющим единичную скорость и направленным под углом атаки а = п/6. Точечный вихрь имеет циркуляцию Г т = 10 и располагается в области течения в точке рт = (2, 1.5)т. На Рис. 1.7, б, показан график точного решения — интенсивности вихревого слоя в зависимости от значений параметра I е [0, 2п).

4 2 0 -2 -4

Y

п 2 п

3 п 2 2 п

-6

аб Рис. 1.7. Линии тока при обтекании эллипса с полуосями а1 : Ь1 =4:1 при наличии вихря в области течения (а) и интенсивность вихревого слоя в зависимости от параметра I (б)

t

На Рис. 1.8, а и б, показаны аналогичные результаты для профиля Жуковского с параметрами а = 3, й = 0.5, Н = 0.3. Скорость набегающего потока, угол атаки и циркуляция вихря те же, что и выше, однако он расположен теперь в точке рт = (1.1, 1)т.

а

3' 0 -3 -6 -9

7*

п 2 п

3 п 2 2 п

б

Рис. 1.8. Линии тока при обтекании профиля Жуковского при наличии вихря в области течения (а) и интенсивность вихревого слоя в зависимости от параметра t (б)

г

Рассмотренная выше методика может быть применена и к рассмотрению так называемых обобщенных профилей Жуковского [75] — крыловых профилей, у которых кромка представляет собой не точку возврата (с величиной угла 0°), а угловую точку с конечным углом. Для этого случая конформное отображение, переводящее внешность круга С радиусом Я на внешность такого крылового профиля К, имеет вид

_2__Л

1 - (Ш)' У

Для задания параметров отображения используются те же соотношения (1.34), что и для классических профилей Жуковского, а также дополнительный параметр а ^ 2. При а = 2 отображение принимает вид классического отображения Жуковского (1.33); в общем случае значение параметра а определяет величину угла на задней кромке в = (2 — а)п.

СШ = 1 аа

Обратное отображение внешности профиля К на внешность круга С радиусом Я имеет более сложный по сравнению с (1.33) вид

С(2) = -Н - а 1 +

2

2 - аа/2\1/а 1 г + аа/2 /

)

и для проведения практических расчетов следует обеспечить выбор непрерывной ветви соответствующей функции.

Обеспечение выполнения условия Чаплыгина — Жуковского в угловой точке приводит к тому же выражению для циркуляции (1.37), что и для классического профиля Жуковского (с учетом отсутствия вихрей в области течения):

Г = -2пЯ|Уто| Бт(а + <р).

(1.38)

Окончательную формулу для вычисления интенсивности вихревого слоя приводить здесь не будем ввиду ее громоздкости, при этом алгоритмический процесс ее получения идентичен вышеописанному. На Рис. 1.9 показаны форма обобщенного профиля Жуковского для параметров а = 3, й = 0.3, Н = 0.25, а =1.9 и распределение интенсивности вихревого слоя 7* на его границе при набегающем потоке единичной скорости, направленном под углом атаки а = п/6, при этом величина 7* стремится к нулю при приближении к угловой точке с обеих сторон.

У

3 0

- 3

- 6 - 9

-12

а б

Рис. 1.9. Линии тока для обобщенного профиля Жуковского (а); интенсивность вихревого слоя в зависимости от параметра I (б)

я 2 Ж

3 я т 2 л

Если задать величину циркуляции отличной от (1.38), положив, к примеру, Г = 0, картина обтекания существенно меняется (Рис. 1.10): угловая точка перестает обтекаться плавно, решение при приближении к ней перестает быть ограниченным, однако остается интегрируемым.

Г

а б

Рис. 1.10. Линии тока при обтекании обобщенного профиля Жуковского при нулевой величине циркуляции (а) и зависимость интенсивности вихревого слоя на профиле от параметра £ (б)

1.5.2. Точные решения трехмерных модельных задач

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

Приведем здесь без подробного описания методики вывода лишь простейшее решение для неподвижной сферы произвольного радиуса Я, расположенной в начале координат, обтекаемой набегающим потоком, имеющим скорость V

В этом случае течение является потенциальным, его потенциал равен

Я3

Ф(г ) = ^ г >(1 + ^3).

Скорость в любой точке области течения, таким образом, может быть вычислена по формуле

v с)=('- цз )v»-(v »■r •

В силу симметрии сферы далее будем без ограничения общности считать, что набегающий поток направлен против оси, проходящей через полюса сферы (от южного с северному), тогда V» = (0, 0, -К»), и введем в пространстве сферическую систему координат, где р — расстояние от центра сферы, р ^ R, 9 — широта, изменяющаяся в диапазоне [0, п] и отсчитываемая от северного полюса к южному, ^ — долгота, принимающая значения в диапазоне [0, 2п). Тогда

r3 I 2рз

Ф(г ) =--К» cos 9,

а компоненты вектора скорости в сферических координатах имеют вид

/ R3 \ / r3 n

Ц1 - p^cos ß, Vv = 0, Vß = 1Ц 1 + 2рз

sin ß.

Тогда интенсивность вихревого слоя на поверхности можно вычислить как y (г) = V (r) х n(r), r е K, в результате получаем

3

Y(r) = 2Я(V» х r),

или, в сферических координатах,

3

Yp = 70 = 0, 7^ = -^К» sin 0.

Поскольку интенсивность вихревого слоя y связана с плотностью потенциала двойного слоя g(r) соотношением (1.6)

— Gradg(r) х n(r) = Y(r),

то можно получить выражение для g(r):

з

g(r) = 3(V^ r),

или, в сферических координатах,

ё = - 2 соб 0.

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

Можно записать также выражение для потенциала скоростей при осесимметричном обтекании эллипсоида вращения, вытянутого вдоль направления набегающего потока, а также «сплюснутого» эллипсоида, у которого фокусы меридионального сечения не лежат вдоль направления набегающего потока [75,93]. Эти задачи решаются, как правило, в эллиптических координатах, и окончательный ответ для интенсивности вихревого слоя имеет весьма громоздкую структуру. Возможно построить и потенциал обтекания произвольного эллипсоида, однако результат в этом случае будет даваться в виде ряда, содержащего специальные функции [62,72]. Отметим также, что при построении решения для потенциального обтекания эллипсоида при решении одной из вспомогательных задач также может быть использован аппарат функций комплексного переменного и конформных отображений [20].

Рассмотрим здесь еще одну задачу, имеющую точное аналитическое решение при наличии завихренности в области течения. Известно, что завихренность в трехмерных задачах должна образовывать замкнутые вихревые линии (вихревые трубки), поэтому наиболее естественно рассматривать замкнутые вихревые кольца. При моделировании нате-кания вихревого кольца на сферу так же, как и в плоском случае при рассмотрении кругового профиля, можно использовать метод отражений, подробно описанный в работе [217] и развитый в работе [288]. Этому вопросу посвящено достаточное число работ, среди которых отметим работы [175,274], однако приведенные там соотношения чрезвычайно громоздки и едва ли годятся для практического использования. Иначе обстоит дело с работой [236], в которой помимо весьма подробного обзора работ по исследованию осесимметричного движения коаксиальных

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

В качестве вихревого кольца рассмотрим тонкую вихревую нить, имеющую в пространстве форму окружности радиусом Я. Будем считать, что ось кольца совпадает с осью Ог, а его центр имеет координату 1, как показано на Рис. 1.11.

Интенсивность вихревого кольца обозначим за Г, это означает, что циркуляция поля скоростей вдоль любого контура, который «протыкается» вихревым кольцом, равна Г.

В соответствии с законом Био — Савара, скорость, индуцируемая вихревым кольцом в точке с радиус-вектором р, равна

Рис. 1.11. Вихревое кольцо

где интеграл берется вдоль вихревого кольца.

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

2ПГ Я2

Ух (2) = Уу (г) = 0, Уг(г) =

(Я2 + (г - г)2)3

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

Уу (г, я, г ) = г( + ^^ > - г),

V уп 1Н2 уп 1Н 2 /

т/ (Г Я г) = Г(2Я2(НХЕХ+Н2Е2) Н\ЕХ -Н\Кх -Н\К2\ Уг(Г ,Я,г) = Г ^ П04 - 4Я2у2 Н'2 Н2' У ,

) (1-39)

где

К'= К(-4Яу). Е,= Е (-Я)

К2 = К( Я ). Е,-Е®

Е2 = Е

12 ' х '12

К(1) и Е(^ — полные эллиптические интегралы 1-го и 2-го рода соответственно [9];

Но = уЯ2+ у2 + (г - г)2, Н\ = \/(Я - у)2 + (г - г)2, Н2 = \/(Я + у)2 + (г - г)2.

Линии тока поля скоростей в плоскости Оуг, индуцируемые вихревым кольцом единичного радиуса (Я = 1), центр которого расположен на расстоянии г = 1 от начала координат и имеющего единичную циркуляцию Г = 1, изображены на Рис. 1.12.

y

HI

-2 -1 ^ ПК

"2

Рис. 1.12. Линии тока поля скоростей, создаваемого вихревым кольцом

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

Z(в)

B2Z

R(B)

b2r

г(в) = _ Г

r2 + z 2

R2 + Z2' R2 + Z2' ~ В

Тогда суммарное поле скоростей, генерируемое двумя кольцами, компоненты которого в плоскости Oyz имеют вид

Vy = Vy(Г, R, Z) + Vy (Г(В), R(B), Z(В)), Vz = Vz(r, R, Z) + Vz(r(B), R(B), Z(B)),

где слагаемые в правых частях равенства выражены формулами (1.39), обладает нулевой нормальной компонентой на поверхности сферы (Vp = 0), а касательная компонента может быть вычислена по формуле

Ve = (Vy cos в - Vz sin в) J z=B cos в.

y=B sin e

Вид линий тока, соответствующих данному полю скоростей, создаваемому рассмотренным выше вихревым кольцом с параметрами R = 1, Z = 1, Г = 1 при наличии обтекаемой сферы с центром в начале координат и радиусом В =1, показан на Рис. 1.13.

z

n

Рис. 1.13. Линии тока поля скоростей, создаваемого вихревым кольцом при его натекании на сферу

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

= Vy C0S ° - V SÍn \z=B cos в.

y=B sin в

Окончательные формулы выписывать здесь не будем ввиду их громоздкости. Зависимость величины от угла в на сфере, отсчитываемого от северного полюса (со стороны которого находится вихревое кольцо) при рассмотренных выше параметрах показана на Рис. 1.14.

о г---—-в

- 3

- 6 - 9

-12

Рис. 1.14. Окружная компонента интенсивности вихревого слоя на сфере в присутствии вихревого кольца

п 4 ^---3 п 4 п

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

1.5.3. Оценка точности приближенного решения

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

Будем считать, что обтекаемый профиль в расчете заменяется N-угольником, вершины которого лежат на границе профиля, а приближенное решение интегрального уравнения представляет собой некоторое распределение, в общем случае кусочно-непрерывное с разрывами в вершинах профиля. Примем, что граница профиля параметризована параметром I, т.е. г = г(¿), I е [0, 2п], а точное решение известно и задано функцией 7*(0.

Наиболее простым и очевидным способом оценки точности численного решения является определение наибольшей по абсолютной величине разности между средними по отрезкам многоугольника величинами интенсивности вихревого слоя 7 и (7*), I = 1, ..., N, вычисленными по приближенному и точному решениям соответственно:

Д7|ся = тах |7 - (7*),! (1.40)

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

Соответствующая ей величина относительной ошибки может быть определена как

max I ^- ( 7 *)«• I

= ^N | ^ . (1.41)

max | (7*)•• | 11

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

Ç 2п

Д7 = | 7(t) - Y*(t) | dlt. (1.42)

Jo

Относительную ошибку приближенного решения можно определить как отношение указанной величины ДY к интегралу от модуля точного решения:

г- 2п

/ | Y(t) - Y*(t)| dlt

¿Y = -. О.43)

/ | Y*(t) | dlt Jo

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

Отметим, что порядок точности численной схемы, рассчитываемый по скорости убывания величины Д7|сл, вычисляемой по формуле (1.40) при к ^ 0, может быть произвольным независимо от распределения интенсивности вихревого слоя по участку границы профиля в точном и приближенном решении. К примеру, кусочно-постоянное численное решение будет обеспечивать в этом смысле нулевую величину ошибки, если оно соответствует средним величинам точного решения. В то же время порядок численной схемы в смысле введенных формулами (1.42) и (1.43) величин абсолютной и относительной ошибок Д7 и в данном контексте существенно более информативен; при этом ясно, что при кусочно-постоянном представлении численного решения порядок может быть не выше первого, при кусочно-линейном — не выше второго. Именно величину £7, определяемую по формуле (1.43), в дальнейшем будем использовать как основную для анализа точности решений.

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

Для этого определим отображение [0, 2п] ^ [0, Ь], где, напомним, отрезок [0, 2п] соответствует диапазону изменения параметра £ при полном обходе исходного контура К; Ь — периметр аппроксимирующего его многоугольника. Введем также естественный параметр р е [0, Ь], величина которого измеряется вдоль границы многоугольника; значения £ = 0 и £ = 2п соответствуют р = 0 и р = Ь.

Обозначив за гI и г¿+1 соответственно радиус-векторы начала и конца ¿-й стороны многоугольника, будем считать, что им соответствуют значения параметра £ = и £ = ¿¿+1 на контуре К и значения параметра р = р1 и р = р1+1 на многоугольнике, I = 1, ..., N. Тогда для любого

значения параметра I е [и, ti+\] можно найти радиус-вектор г(¿) точки на исходном контуре и значение точного решения для интенсивности вихревого слоя. Чтобы найти численное решение, соответствующее точке г^), опустим перпендикуляр из этой точки на сторону многоугольника (Рис. 1.15) и найдем значение параметра рследующим образом:

т= + (г(0 - г^ ■ (г,+, - гi). (1.44)

|Гi+1 - Гi 1

г 1+\ = г(и+\)

Р(+\

Р

Рис. 1.15. Проекция точки на границе исходного профиля на сторону многоугольника

В результате можно считать, что зависимость 7(0 = 7(р(0), где 7(р) — распределение интенсивности вихревого слоя, найденное в расчете на границе многоугольника, а р(0 определяется формулой (1.44), определяет правило проецирования численного решения на исходный контур К. Это позволяет вычислять погрешность по формуле (1.42), используя в ней в качестве численного решения 7(0 зависимость 7(0.

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

Для более точного описания формы границы обтекаемого профиля ее можно заменить не многоугольником, сторонами которого являются прямолинейные отрезки, а фигурой более сложной формы, которую назовем криволинейным многоугольником. Пример подобного алгоритма аппроксимации представлен в работах [68,83,188,195,225]. В этом случае все вышеприведенные рассуждения остаются в силе за исключением сравнительно простого соотношения (1.44); вместо него для вычисления проекции численного решения, заданного на криволинейном многоугольнике, на исходный контур также требуется решение нелинейного уравнения.

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

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

1.6. Выводы по главе 1

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

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

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

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

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

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

Глава 2. Расчетные схемы для решения плоских

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

Рассмотрим способы построения численных схем решения граничных интегральных уравнений, соответствующих N - и Т-моделям в плоских задачах моделирования обтекания профилей, рассмотренным в разделах 1.2 и 1.4. Соответствующие схемы назовем N - и Т-схемами, при этом под «схемой» будем понимать совокупность методов численного решения соответствующих задач и реализующих их вычислительных алгоритмов.

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

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

Для анализа точности существующих и разрабатываемых схем будут использованы точные решения модельных задач, построенные в разделе 1.5.

2.1. Классическая схема метода дискретных вихрей

Рассмотрим классическую N-модель, выражаемую интегральным уравнением (1.11) (в двумерном случае при d = 2) и совпадающим с ним интегральным уравнением (1.25), которое можно записать в виде

£KPn(r, 0)y(0) dit = fn(r), r G K, (2.1)

где

Pn(r, 0 = -T(r) • Q(r - 0 = -T2rn)|r (-00, (2.2)

t(r) — орт касательной в соответствующей точке границы профиля, выбираемый таким образом, что при обходе контура K в его направлении область течения F остается справа; Q (r - 0) — взятый с обратным знаком градиент фундаментального решения уравнения Лапласа (в двумерном случае).

Для простоты в этом разделе будем полагать, что обтекаемый профиль неподвижен (UK(r) = 0, r G K), тогда правая часть интегрального уравнения (2.1) принимает более простой по сравнению с (1.24) вид:

fn(r) = f (r) • n(r) = -n(r) • (V+ jf (k x Q(r - 0)0(0 dS^j =

= - (n(r) • Vœ + jf Pn(r, 0)0(0) dS^,

где n(r ) — орт внешней по отношению к профилю нормали в точке r ; Vœ — скорость набегающего потока; k — орт нормали к плоскости течения; 0(0 — распределение завихренности в области течения F.

Завихренность, как было отмечено в разделе 1.3, в плоских задачах естественно считать скалярным полем, поскольку вектор завихренности имеет лишь одну ненулевую компоненту в направлении нормали к области течения. Распределение завихренности в области течения 0(p) будем считать заданным при помощи множества дискретных вихревых элементов с известными циркуляциями Г w и положениями pw, w = 1, ..., N0:

No

0(p) = У Г wP - Pw) ,

w=1

где £(p) — двумерная дельта-функция Дирака [24].

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

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

N

7(г) = Е - Р)), г е K, (2.3)

)=1

где Pj — точки расположения вихревых элементов, задаваемые при построении расчетной схемы по одной на каждую панель; г|7) — циркуляции вихревых элементов, подлежащие определению (верхний индекс (7) введен, чтобы подчеркнуть их связь с вихревым слоем и одновременно отличать данные вихри от тех, что моделируют распределение завихренности в области течения).

С учетом представления интенсивности вихревого слоя (2.3) уравнение (2.1) можно записать в виде

N

Е Pn(r, Р))г(7)= Ы(Г), Г е K. (2.4)

)=1

Далее на каждой панели профиля выбирают по одной так называемой контрольной точке 01, I = 1,..., N, — точке коллокации, в которой обеспечивается выполнение равенства (2.4). В результате получается система

линейных алгебраических уравнений вида

N

Е РпО, Р) )ГУ(7)= Ш), I =1,..., N, (2.5)

)=1

где

( ^ А

^(01) = - ( а(01) • V+ Е Рп(01, Рт)Гда ] .

^ т=1 '

Величина Рп(г, 0), определяемая формулой (2.2), выражает собой нормальную компоненту скорости в точке г, индуцируемой вихревым элементом единичной циркуляции, находящимся в точке 0 Ясно, что Рп(г, 0) ^ то при |г — 0| ^ 0, поэтому для ограничения скоростей, индуцируемых вихревыми элементами, вводят понятие радиуса вихря, называемого также радиусом дискретности, — малую положительную величину е с размерностью длины и вместо функции Рп(г, 0), заданной формулой (2.2), рассматривают функцию

Ртсг, о = — 2 т(г' ;<г — ? 21. (2.6)

2пшах{|г — о|2, е21

Отметим, что регуляризация скорости, индуцируемой вихревым элементом и выполняемая по правилу (2.6), соответствует переходу от модели точечного вихря (вихревой нити) к более сложной модели вихря.

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

Если же рассмотреть круглый вихрь радиусом е, содержащий ту же суммарную завихренность , которая теперь распределена равномерно по его площади (такой вихрь обычно называют круглым вихрем Рэнки-на [72,92,93]), то индуцируемая им скорость имеет вид

Vе (р) = ГК(р) = Г--** Р п . (2.7)

2пшах{ |р|2, е2)

Сопоставление выражения (2.7) с (2.6) с учетом вышеуказанного смысла величины Рп(г, 0) показывает, что замена функции влияния Рп(г, 0) на Рп(г, 0) соответствует переходу от представления завихренности точечными вихрями к ее представлению круглыми вихрями Рэнкина.

Отметим, что описанная методика «сглаживания» особенности является наиболее простой, но не единственно возможной. В монографии [147] приведено не менее четырех способов регуляризации поля скоростей с подробным обсуждением влияния выбора функции сглаживания на точность результатов; этому же вопросу посвящена работа [139] и большое количество других публикаций, главным образом, зарубежных; обзор ранних исследований в этом направлении сделан в [264]. При этом вопрос выбора величины радиуса сглаживания е также является дискуссионным.

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

Рис. 2.1. Схемы расстановки точки расположения вихрей и контрольных точек в методе дискретных вихрей

При этом система уравнений (2.5) оказывается близкой к вырожденной (как обсуждалось в разделе 1.3, уравнение (2.1) имеет бесконечное множество решений), поэтому для выделения единственного решения требуется дополнительное условие, которое получается из (1.27) подстановкой туда представления решения в виде (2.3):