CC BY
43
УДК 532.5
DOI: 10.18698/1812-3368-2016-6-93-109
ЧИСЛЕННАЯ СХЕМА ВЫСОКОГО ПОРЯДКА ТОЧНОСТИ ДЛЯ ОПРЕДЕЛЕНИЯ ИНТЕНСИВНОСТИ ВИХРЕВОГО СЛОЯ ПРИ РЕШЕНИИ ДВУМЕРНЫХ ЗАДАЧ АЭРОГИДРОДИНАМИКИ ВИХРЕВЫМИ МЕТОДАМИ
К.С. Кузьмина И.К. Марчевский
kuz-ksen-serg@yandex.ru iliamarchevsky@mail.ru
МГТУ им. Н.Э. Баумана, Москва, Российская Федерация
Аннотация
Рассмотрены вопросы развития вихревых методов для численного моделирования двумерного обтекания профилей потоком вязкой несжимаемой среды. Используемые в вихревых методах расчетные схемы и алгоритмы определения интенсивности вихревого слоя предполагают аппроксимацию профиля многоугольником, состоящим из прямолинейных отрезков-панелей, интенсивность вихревого слоя на которых принимается кусочно-постоянной. При применении оптимальных схем данный подход позволяет обеспечивать погрешность от O(h2) до O(h3) в зависимости от рассматриваемого профиля (h — длина панели). Разработана новая расчетная схема, в которой интенсивность вихревого слоя полагают кусочно-линейной, либо кусочно-квадратичной функцией, а также учтена криволинейность поверхности профиля. Для получения основной системы линейных алгебраических уравнений использован метод наименьших квадратов вместо применяемого условия типа коллокаций в контрольных точках или в среднем на панелях. Показано, что разработанная схема имеет более высокую точность по сравнению с точностью ранее известных схем: для рассмотренных тестовых задач (обтекание круга, эллипса, профиля Жуковского) схема обеспечивает погрешность на уровне O(h5)
Ключевые слова
Вихревой метод, вихревой слой, схема высокого порядка, метод наименьших квадратов, кривизна, кусочно-полиномиальная аппроксимация
Поступила в редакцию 25.05.2016 © МГТУ им. Н.Э. Баумана, 2016
Работа выполнена при финансовой поддержке гранта Президента России для молодых российских ученых — кандидатов наук (проект MK-7431.2016.8)
Введение. Вихревые методы, относящиеся к классу бессеточных лагранжевых методов вычислительной гидродинамики, — весьма мощный инструмент во многих инженерных приложениях, в частности, при решении сопряженных задач гидроупругости, когда течение среды можно полагать несжимаемым, а форма расчетной области изменяется вследствие движения помещенных в нее тел,
обусловленного, в свою очередь, действующими на них гидродинамическими нагрузками.
При рассмотрении двумерных задач известно несколько подходов к решению уравнений Навье — Стокса вихревыми методами [1-5]; наиболее удобным для практического использования является метод вязких вихревых доменов (ВВД) [5], в котором для моделирования эволюции завихренности в вязкой жидкости используется так называемое поле диффузионных скоростей [3, 6]. Точность моделирования течения и расчета гидродинамических нагрузок определяется множеством факторов, наиболее существенные из которых следующие:
1) точность аппроксимации профиля;
2) точность определения интенсивности вихревого слоя на поверхности профиля;
3) точность аппроксимации вихревого следа и моделирования его эволюции.
В рамках настоящей работы рассмотрим возможные пути повышения точности за счет первых двух факторов; они представляются наиболее значимыми, поскольку источником формирования вихревого следа является именно вихревой слой, генерируемый на поверхности профиля.
Как правило, в реализациях вихревых методов, применяемых при решении двумерных задач гидродинамики и гидроупругости, интенсивность вихревого слоя на профиле находят как решение сингулярного граничного интегрального уравнения первого рода [1]. Опыт показывает, что такой подход может приводить к существенным погрешностям, а в ряде случаев — к получению качественно неверного решения [7]. Это особенно актуально при расчете обтекания гладких профилей, однако может быть важно и при расчете течений вблизи острой кромки [8].
Известен также альтернативный подход, предполагающий сведение задачи о поиске интенсивности вихревого слоя на профиле к решению интегрального уравнения типа Фредгольма второго рода [9]. Эффективность применения этого подхода для решения практических задач рассмотрена в работах [7, 8, 10, 11], там же описаны расчетные схемы, основанные на использовании указанного подхода, и показана возможность существенного повышения точности решения задач.
Дальнейшее повышение точности решения ограничено качеством аппроксимации обтекаемого профиля. В настоящей работе предложен алгоритм и построены необходимые квадратурные формулы, позволяющие развить идеи подхода [9] с учетом кривизны панелей профиля. В рамках предложенного подхода искомая интенсивность вихревого слоя может быть кусочно-постоянной, кусочно-линейной, либо кусочно-квадратичной функцией на каждой панели (в рамках традиционного подхода интенсивность вихревого слоя — кусочно-постоянная функция на прямолинейных панелях, аппроксимирующих профиль).
Отмеченные пути повышения точности моделирования могут показаться аналогичными тем подходам, которые применяют в так называемых панельных
методах [12], однако между ними есть существенные различия. Прежде всего, в настоящей работе решение не предполагается непрерывным (и тем более непрерывно дифференцируемым) вдоль профиля — это необходимо для корректного моделирования обтекания профилей с угловыми точками и острыми кромками. В связи с этим для вычисления интегралов вдоль криволинейных участков профиля используют формулы численного интегрирования Гаусса вместо интегрирования разложений в степенные ряды. Кроме того, для получения основной системы линейных алгебраических уравнений использованы идеи метода наименьших квадратов: истинное решение должно минимизировать интеграл, вычисляемый по твердым границам области течения («стенкам»), от квадрата разности скорости среды, определяемой искомым распределением завихренности, и заданной скорости обтекаемого профиля.
Описанный подход и построенный на его основе численный алгоритм позволяют значительно повысить точность определения интенсивности вихревого слоя в вихревых методах.
Постановка задачи. Течения вязкой несжимаемой среды описывают уравнениями Навье — Стокса
- дV - - - Vp
V-V = 0, -+ (У-У)У = уДУ--¡-,
дt р
где У(г, t) — поле скоростей; V — коэффициент кинематической вязкости; р — плотность среды, принимаемая постоянной; р(г, t) — давление.
Рассмотрим задачу о внешнем обтекании профиля неограниченным потоком среды, тогда граничными условиями являются условие прилипания жидкости на профиле и условия затухания возмущений на бесконечности:
У(-, 0 = Ук(-, t), - е К; У(-, t) ^ У», р(-, t) ^ р», | г|^>а>.
Здесь Ук (г, t) — скорость точек профиля, которая предполагается известной.
Уравнения Навье — Стокса могут быть записаны в форме Гельмгольца с использованием поля завихренности О(г, t) = Ух У (г, t):
дО. - -
-+Ух (Охи ) = 0, (1)
дt
где и (г, I ) = У (г, I) + Щ (г, I), Щ (г, t) — так называемая диффузионная скорость, пропорциональная коэффициенту вязкости среды и определяющая эволюцию завихренности в вязкой жидкости [3],
т-,- ч (УхО) хО
Щ (г, г ) = V--^-.
IО |2
Согласно уравнению (1), завихренность, имеющаяся в начальный момент времени в области течения, движется со скоростью и, в то время как «новая»
завихренность генерируется только на обтекаемом профиле. Примем, что начальное распределение завихренности в области течения Q(r, t0) известно.
Влияние обтекаемого профиля на течение эквивалентно суперпозиции влияний присоединенного вихревого слоя с интенсивностью yatt (r, t), присоединенного слоя источников интенсивностью qatt (r, t) и свободного вихревого слоя неизвестной интенсивности y(r, t). Эти слои располагаются на обтекаемом профиле, присоединенные слои моделируют движение и, в том числе, возможное деформирование профиля, поэтому их интенсивность определяется скоростью движения точек профиля:
Уatt (r, t) = Vk (r, t) -x(r, t), qatt (r, t) = Vk (r, t) • n(r, t), r e K,
где n(r, t), x(r, t) — орты нормали и касательной к профилю [5, 7].
В настоящей работе для простоты рассмотрим обтекание неподвижного не-деформируемого профиля, поэтому yatt (r, t) = 0, qatt (r, t) = 0, однако это предположение не является принципиальным и может быть опущено. Примем также, что скорость набегающего потока является постоянной V» = const.
По известному распределению завихренности с помощью закона Био — Савара может быть восстановлено поле скоростей среды
лтп а лт ^ 1 r^(?, *) х (r - ^ ^ f(?, t) х (г - ?)
V (r, t ) = У»+—J-———-dS+—ф-———-dls. (2)
2% S | r - s |2 2% K | r - s |2
Здесь S — область течения; K — обтекаемый профиль; у = yk, Q = Qk — векторы интенсивности вихревого слоя и завихренности в области течения; k — орт нормали к плоскости течения, П(г ) х x(r) = k.
Интенсивность вихревого слоя y(s, t) можно определить из граничного условия на неподвижном профиле V(r, t) = 0, r e K.
Рассмотрим простейшую модельную задачу и примем, что завихренность в области течения отсутствует (Q(r, t ) = 0), и в этих условиях требуется определить интенсивность вихревого слоя на профиле. Математически указанная задача эквивалентна расчету обтекания профиля потоком идеальной несжимаемой среды при заданной циркуляции вокруг профиля. При расчете течений вязкой несжимаемой среды свободный вихревой слой сходит в поток на каждом временном шаге, формируя вихревой след вблизи и позади профиля, а описанную выше задачу расчета интенсивности генерируемого на поверхности вихревого слоя решают на каждом шаге расчета.
Интегральное уравнение для определения интенсивности вихревого слоя. Учитывая, что неизвестная интенсивность y(s, t) соответствует свободной завихренности, являющейся частью вихревого следа в области течения, в соответствии с (2) можно показать, что предельное значение скорости среды со
стороны профиля имеет вид (здесь и далее зависимости всех величин от времени опущены)
т7/-\ т7 1* х(---)лт 1^(-) -/-О - ^
У- (г) = У» +—С—---— --х п(г) , г е К. (3)
2л к | г - 5|2 ^ 2 )
«Классический» подход, обычно используемый в реализациях вихревых методов, предполагает, что неизвестную функцию у(г) следует определять из условия равенства нулю нормальной компоненты предельного значения поля скоростей на профиле:
У-(-)-П(-) = 0, - е К. (4)
С учетом (3) относительно искомого распределения у(г) получается интегральное уравнение первого рода, которое является сингулярным; интеграл в нем следует понимать в смысле главного значения по Коши [1]. Численное решение этого уравнения может приводить к значительным погрешностям, а в ряде случаев — и к получению качественно неверного результата.
В то же время известен альтернативный подход, из которого, в частности, следует [9], что условие (4) эквивалентно условию
У-(-)-т(-) = 0, - еК, (5)
которое означает равенство нулю касательной компоненты предельного значения поля скоростей на профиле.
Можно показать, что (5) приводит к интегральному уравнению второго рода типа Фредгольма
—С ,- -,2— Т(5)^--= -У»-х(г). (6)
2л к | г - 512 2
Ядро полученного уравнения является ограниченной функцией для случая гладкого профиля:
| П(г ) - (г - 5 )| к (г)
Цш —— -2— ,
|---|^0 | г - 512 2
где к (г) — кривизна профиля в соответствующей точке.
Интегральное уравнение (6), как и сингулярное уравнение, следующее из (4), имеет бесконечное множество решений. Для выделения единственного решения их требуется решать совместно с дополнительным условием, задающим величину суммарной циркуляции Г скорости вокруг профиля, которая обычно известна из постановки задачи:
|у(-№ = Г. (7)
к
Для профилей простейших форм можно применить известный метод конформных отображений и решить задачу определения интенсивности вихревого слоя аналитически. Методика получения точных решений для кругового,
эллиптического профиля и профиля Жуковского описана в работе [13]. Эти решения будут использованы в качестве «эталонных» для оценки точности построенных расчетных схем.
Расчетные схемы для определения интенсивности вихревого слоя. Рассмотрим две расчетные схемы для определения интенсивности вихревого слоя.
Схема с прямолинейными панелями. Как было отмечено выше, в вихревых методах обтекаемый профиль обычно аппроксимируют ^-угольником, стороны которого имеют длины Ь, и называются «панелями», а интенсивность вихревого слоя постоянна на панелях.
В рамках сделанных предположений интеграл в уравнении (6) можно заменить суммой интегралов по отдельным панелям, которые будут пропорциональны соответствующим интенсивностям. Точность решения задачи можно значительно повысить, если потребовать выполнения уравнения (6) не в отдельных точках на профиле, называемых обычно «контрольными», а в среднем на панелях
1 £у, jij ^
2nLi
■i j—1 Ki
V Kj
I r - s I
dlr — Уда
2
i — 1, ...,N.
Коэффициенты выписанной системы линейных алгебраических уравнений (СЛАУ) могут быть вычислены аналитически [10, 11].
Вычислительные эксперименты показывают, что предложенная схема обеспечивает для различных профилей погрешность определения циркуляций вихревых элементов, сходящих с профиля в поток, порядка 0(Н2) или 0(Н3), где Н — средняя длина панели (рис. 1).
100 200 а
Рис. 1. Погрешность определения циркуляций вихревых элементов (в логарифмическом масштабе) для схемы с прямолинейными панелями: а — эллиптический профиль (штриховая линия соответствует третьему порядку точности); б — профиль Жуковского (штриховые линии соответствуют второму и третьему порядку
точности)
Схема, учитывающая криволинейность профиля. При разработке расчетной схемы повышенной точности, во-первых, следует учесть криволинейность профиля, а во-вторых, повысить точность аппроксимации интенсивности
вихревого слоя на каждой панели, сделав соответствующее распределение кусочно-линейным или кусочно-квадратичным.
Учет каждого фактора в отдельности позволяет обеспечить лишь незначительное повышение точности и не повышает порядок точности схемы.
Известно [1], что интенсивность вихревого слоя неограниченно возрастает вблизи «внешних» угловых точек профиля (т. е. в тех случаях, когда угол между касательными к сторонам угла больше развернутого, если смотреть со стороны потока), поэтому гладкие части профиля при построении схем повышенного порядка точности с необходимостью следует аппроксимировать гладкими кривыми. Будем предполагать, что форма профиля известна точно и она описывается кусочно-гладкими параметрическими зависимостями х = х(£), у = у($), I £ [0,2 л). Таким образом, можно считать, что известны не только координаты отдельных точек на профиле, определяющих концы панелей, но и направления касательных к профилю в этих точках.
Потребуем, чтобы кривая, аппроксимирующая профиль, проходила через заданные точки, а также имела в этих точках то же направление касательной, что и исходный профиль.
Обозначим начало и конец г-й панели через Сг и Сг+1 (эти точки соответствуют значениям параметра I = ti и t = в параметрических уравнениях
профиля); пусть хг° — единичный вектор, сонаправленный с вектором СС+ъ и® — единичный вектор, ортогональный вектору хг° (рис. 2). Под «панелью профиля» теперь будем понимать прямолинейный отрезок, соединяющий точки С и С+1; при этом в общем случае участки профиля имеют с соответствующими панелями лишь общие концы.
На каждой панели введем локальную систему координат Ci ^ ^, тогда точки Ci и Ci+i имеют локальные координаты (0,0) и (Li ,0), где Li =| QQ+i|. Для построения интерполяционной кривой сначала вычислим тангенсы углов наклона касательной к панели и ^ (см. рис. 2), используя следующие соотношения:
/
О х
Рис. 2. Кусочно-полиномиальная аппроксимация профиля
x+ (ti )sin9,' - y + (ti )cos 9,' x- (ti+1)sin9,' - y- (ti+1)cos 9; tg Ф* — —TTTs-~ , . n ; tg Vi — ■---
X + (^ )С08 9; + У + (^ )в1И 9; X(^ + 1)С08 9; + у- (^ +1)81и 9;
Здесь 9; — угол между г-й панелью и осью Ох; х+ , х-, у+, у'- — право- и левосторонние производные параметрических зависимостей х($) и у($), вычисляемые в соответствующих точках; фг и у, могут принимать значения из промежутка (-71 /2; л /2).
Уравнение, задающее положение интерполяционной кривой, аппроксимирующей исходный профиль, в локальной системе координат над ;-й панелью имеет вид
р ф = ¿(¿Ч) г ^ +
Ь; V Ь;
где условия рг (0) = 0, рг (Ь) = 0 выполняются автоматически, коэффициенты а; и Ь; определяют из условий р-(0) = гдф;, р'(Ьг-) = ^у: а, ^ф;, Ь = -гдф;.
С учетом изложенного положение произвольной точки М, лежащей на интерполяционной кривой и имеющей координату (абсциссу) в локальной системе координат, задает радиус-вектор ОМ (¿) = 0С; + р; (¿)й°.
Если участок исходного профиля над 2-й панелью представляет собой гладкую кривую класса С4, погрешность аппроксимации профиля будет величиной порядка 0(Ь4). Этот результат легко получить, построив соответствующие разложения по формуле Тейлора.
Аппроксимируем неизвестное распределение интенсивности вихревого слоя на участке профиля над 1-й панелью квадратичной функцией
У; (¿) = +Р; ^ + 5; (8)
Ь Ь
Неизвестные коэффициенты а;, Р;, 5;, I = 1,..., Ы, определяют из интегрального уравнения (6) с дополнительным условием (7).
Для отыскания приближенного решения уравнения (8) при принятой зависимости у; (¿) используем метод наименьших квадратов и поставим задачу безусловной минимизации функции
* — if-L | Щ-lll у(г Ж -lf>+V.-X(r)'
J 2l K I r - s I2 2
K
v
dlr -
-xf I y(r)dlrmin (9)
по всем возможным значениям параметров ai, ßi, , X.
И внутренний, и внешний интегралы в (9) можно заменить суммами интегралов вдоль криволинейных участков профиля над панелями:
N
(
Ш V Г 1 N г n(r) •(r ~5) п\А1 У(г) + т? •
Ш = XI — X J -„— y(s )dls —— + V® • т(г)
г=1 к [ 2л j=i к, I r ~ 5I2 2
V
dL ~
f
N
Л
Xjy(r) dlr ~Г
Vi=1 к
► min.
(10)
Все интегралы вычисляют в локальных системах координат, при этом для участков профиля над г-й и ;'-й панелями имеем г = Г (5) и 5 = 5; (С):
Г (5) = ос, + 5х,0 + р (5)йг0, щ (С) = ос;+Сх° + р (дя°.
Якобиан преобразования координат при переходе к локальным координатам имеет вид
ь (5) = ^т 1м (5)=>Д+(Р^))2,
и в результате, обозначая для простоты у(Г (5)) = у г (5), я(Г (5)) = Я,- (5), х(Г (5)) = х, (5), минимизируемую функцию (10) можно записать в виде
N Lj( ' " Lj
Ш = X I — X J n®^®~(0)
i=1 0 4
2л 0 I r© ~ 5,(C)|2
i N Li
у, (QJj (Qdi-^+V®^)
Ji
Л
X I yi (^)Ji Й№-Г
Vi=1 0
min.
(11)
Поскольку информация об исходном профиле, согласно сделанным предположениям, ограничивается координатами узлов и направлениями касательных в них, орт вектора касательной х, (5) в формуле (11) вычисляют не для исходной кривой, а для интерполянта
- 0 i/s-ч-0
Ti + pi&n
I-0 i/s-ч-0
I Ti + pi^H
Погрешность аппроксимации орта касательной составляет величину порядка 0(1^) для кривых класса С4. Орт вектора нормали Я (5) выбирают ортогональным к х, (5).
Подставляя в (11) выражение для искомого распределения интенсивности вихревого слоя (8) и вводя обозначения
| Я (5)'(Г(5) ~ (С)) М = I;) (5), ) л (5)5^5 = V,
получаем функцию, зависящую только от неизвестных коэффициентов а к, Рк, 5к, к = 1, ..., N, и множителя Лагранжа X:
n цг ( 1 n
Щ = — Е(() + Р^(Ц) + 3 ^(Ц))
i =1 0 ^ ¿ = 1
- 2 (а +рг-^+5г Ц- 1 +V» • хг (Ц) 1 т (Ц) йЦ-
/ n
Ы0) + J + J )-Г
V i=1
■ min.
Минимальное значение этой функции достигается в случае равенства нулю частных производных по всем переменным
ящ N Цг ( 1 N ,
^ = Е |2 Де (а ()+Р ¿Ц)+ 5 ¿1^))-
аак г=10 ( 2л ¿=1
-1 (аг + Рг + 8г Ц-1 + V» • Хг (Ц) 1 (^^ -1I (Ц)^ - Т0 = 0; 2 ( Ц Ц) Д 2л 2 I
ЯЩ N ¿г А 1 n
= S J2 fE(а,7f(!) + Р,7)?>(!) + 5<»(!))-
а г+рД + 8^ | + V^-t,(!) 2 V L L2
( TL® -1 V 2л 2 Lk у
Ji (!)d!-Ui1) =0;
ЯЩ N Li ГА 1 n ,
|Щ = S12 II fs (а ^^ f(!) + Р 11)(!) + 5 jl «(!))-
a8k г=1 0 VV 2л Д - -1 Ааг +рг-! + 8г ) + ^ (!)
V 2л 2 L2k J
Ji (!)d! - X/k2) = 0;
ЯЩ = Е ( I (0) + Рг I (1) +8г I (2) )-Г = 0.
ЯХ г=1
Полученные выражения могут быть записаны в более простой форме. Если обозначить интегралы от известных функций через
J(ii= 1 iip}(!)i{s(!)Jm (!)d!;
0
1тп J -^тп(Ц) Тт
0 Цт
неизвестные величины — через у(м), м = 0,1,2, где у(0) = а;, у1;1'= Р;,
у® = 5;, = 1,..., N, то записанная выше СЛАУ может быть представлена
в виде
N 2
XXyf> tVXJT X J(r^ JT +TJ("+r^ ~Jcr) =
^ 1 N ( ^ 1 N
4 j=1u=0
( ir \ 1 1 ir\
V 2л2 i=1 j 2л i=1
T(u,r
j T1
V Lk У
r
1 -(r,u) , 1 t(u+r) ^
2л jk 2Jj Lk
Ar) _
NLi . _ . ^ 1 5 Л
= -Х|(-хг(5))|-(5)-5- ¿5, ^ = 1, •••,N. Г = 0,1,2,
,=10 V- Чк
N 2
хху;м) /(и) =г. (12)
;=1и=0
Если рассматривать кусочно-линейное распределение интенсивности вихревого слоя на участках профиля над панелями, то переменные Г и и (индексы, по которым выполняется суммирование) в системе (12) будут принимать значения 0 и 1.
Для приближенного вычисления всех интегралов, входящих в систему (12), были применены квадратурные формулы Гаусса наивысшего порядка точности, использующие Пр узлов, в соответствии с которыми
w b~аП^Р Ja + b b~a „ ^ J f Шх *——X®kf I——+——Ck I,
a 2 k=1 V 2 2 У
где значения весовых коэффициентов ю и положения гауссовых точек Хк, к = 1,..., Пр, выбирают стандартным образом [14].
В настоящей работе использовано значение п^ = 7. В тех случаях, когда приближенные формулы Гаусса не обеспечивали достаточной точности вычисления интегралов, что обычно наблюдалось при вычислении интегралов по соседним панелям, применяли процедуру дополнительной разбивки каждой панели на Nadd = 4 панели меньшей длины. При необходимости вновь образованные «мелкие» панели подвергали повторной доразбивке.
Вычислительные эксперименты показывают, что учет якобиана /^ (5) мало влияет на точность решения при достаточно большом числе панелей, поэтому на практике можно приближенно принять, что /^ (5) = 1; это несколько упрощает расчеты.
Отметим, что при реализации описанного алгоритма в каком-либо математическом пакете (МаНАВ, МаШешаИса и др.) могут быть использованы встроенные алгоритмы численного интегрирования, что может негативно влиять на время вычисления квадратур, однако значительно упрощает разработку соответствующих подпрограмм.
Вычислительный эксперимент. Рассмотрим задачу расчета обтекания кругового профиля радиусом Я = 1, эллиптического профиля с полуосями а1 = 1,0, Ь1 = 0,5 и профиля Жуковского с параметрами а = 3,5, а = 0,4, к = 0,3, установленного под углом атаки р = - /6.
Точное аналитическое решение. Точное аналитическое решение задачи для эллиптических профилей и профилей Жуковского, используемое для контроля
точности, находят методом конформных отображений [13, 15], и для интенсивности вихревого слоя на профиле можно записать формулу
у. й )=2У„ 5щ(ф+р->) + г/№), н = А (13)
1- «2
№ а _ф) + н )2
Здесь I е [0; 2л) — параметр, определяющий положение точки на профиле; параметрические зависимости, задающие форму профиля, имеют вид х(£ } = Ие 2 (£), у(£ ) = 1т 2^), где
г«) = х() + iy(t) = 1 )|, Х(0 = Яе'^_ф) + Н. (14)
21 ))
Для эллиптического профиля следует принять
а = л]а2 _ Ь2, Я = я1 + Ь1, ф = 0, к = й = 0,
где а1, Ь1 — большая и малая полуоси эллипса; для профиля Жуковского к
Я =| Н _ а |, ф = агС^—, а параметры а, й и к определяют длину, толщину и а
кривизну профиля.
Для эллиптического профиля циркуляцию скорости Г формально можно выбрать произвольно, в рассматриваемой модельной задаче примем Г = 0 независимо от угла атаки; для профиля Жуковского циркуляцию скорости Г выбираем из условия конечности скорости потока на кромке: ее значение пропорционально скорости набегающего потока и зависит от формы профиля и его угла
атаки: Г = _2лУт 8т(р + фХл/к7 + а2 + й).
Численное решение и оценка его точности. С использованием разработанной расчетной схемы, учитывающей криволинейность участков профиля над панелями, сначала вычисляли распределение интенсивности вихревого слоя по панелям (находили значения коэффициентов а,, р, и 5, ), а затем выполняли его интегрирование вдоль панелей, чтобы получить суммарную циркуляцию вихревого слоя на каждой панели
Г, = |Га,- + р, Ц- + 5,-У],($)й\, I = 1,...,N.
0 V Ц Ц)
Аналогичную величину Г. определяли путем интегрирования точного решения (13):
^+1 ,-
г. = I у*(0ч/х'(02 + у'(02йл, I = 1,..., N
ч
где значения параметра ti и ti+1 соответствуют началу и концу г-й панели, зависимости х^) и у(t) находят из формулы (14).
Погрешность численного решения определяется как наибольшее по модулю отклонение между точным и численным решением, рассчитанное по суммарной
циркуляции вихревого слоя на панелях профиля: ДГ = тах |г, - Г* |.
г
Результаты расчетов показывают примерно одинаковую погрешность для рассмотренных модельных задач в предположениях о кусочно-линейном и кусочно-квадратичном распределении интенсивности вихревого слоя по участкам профиля. Следует отметить, что в случае кусочно-квадратичной аппроксимации матрица системы (12) становится плохообусловленной.
Во всех рассмотренных примерах (для кругового, эллиптического профиля и профиля Жуковского) погрешность численного решения в указанном выше смысле имеет порядок 0(Н5) (рис. 3).
ДГ
ю-3 10^ 10"5 10"6 10"7 ю-8
ДГ
10"1 10"2 10"3 10"4
10 20 40 80 N в
Рис. 3. Погрешность определения суммарной циркуляции на участках профиля над панелями (в логарифмическом масштабе) при использовании разработанной
расчетной схемы (12): а — круговой профиль; б — эллиптический профиль; в — профиль Жуковского (штриховая линия соответствует пятому порядку точности)
Заключение. Для моделирования течений вязкой несжимаемой жидкости вихревыми методами разработана новая расчетная схема повышенной точности для определения интенсивности вихревого слоя на обтекаемом профиле. Повышение точности решения соответствующего интегрального уравнения обеспечивается корректным учетом криволинейной формы профиля. Распределение завихренности принято кусочно-линейным, либо кусочно-квадратичным, при этом схема допускает разрывы решения на границах панелей, что позволяет
корректно моделировать обтекание профилей с острыми кромками, в частности, профиля Жуковского. Вычислительные эксперименты на модельных задачах, имеющих точное аналитическое решение, показывают, что предложенная схема обеспечивает пятый порядок точности.
ЛИТЕРАТУРА
1. Lifanov I.K., Belotserkovsky S.M. Methods of discrete vortices. CRC, 1992.
2. Cottet G.-H., Koumoutsakos P. Vortex methods: theory and practice. Cambridge University Press, 2000.
3. Дынникова Г.Я. Движение вихрей в двумерных течениях вязкой жидкости // Известия РАН. Механика жидкости и газа. 2003. № 5. С. 11-19.
4. Калугин В.Т., Мордвинцев Г.Г., Попов В.М. Моделирование процессов обтекания и управления аэродинамическими характеристиками летательных аппаратов. М.: Изд-во МГТУ им. Н.Э. Баумана, 2011. 527 с.
5. Андронов П.Р., Гувернюк С.В., Дынникова Г.Я. Вихревые методы расчета нестационарных гидродинамических нагрузок. М.: Изд-во МГУ им. М.В. Ломоносова, 2006. 184 с.
6. Ogami Y., Akamatsu T. Viscous flow simulation using the discrete vortex model — the diffusion velocity method // Comput. and Fluids. 1991. Vol. 19. Nc. 3/4. P. 433-441.
7. Kuzmina K.S., Marchevsky I.K. On numerical schemes in 2D vortex element method for flow simulation around moving and deformable airfoils // Proceedings of Summer School-Conference "Advanced Problems in Mechanics 2014". St. Petersburg, 2014. P. 335-344. URL: http://www.ipme.ru/ipme/conf/APM2014/2014-PDF/2014-335.pdf
8. Кузьмина К.С., Марчевский И.К., Морева В.С. О точности расчетных схем вихревых методов при моделировании обтекания профилей с угловой точкой // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. № 2. C. 234-249.
DOI: 10.7463/0215.0756954 URL: http://technomag.neicon.ru/doc/756954.html
9. Kempka S.N., Glass M.W., Peery J.S., Strickland J.H. Accuracy considerations for implementing velocity boundary conditions in vorticity formulations. SANDIA REPORT SAND96-0583 UC-700, 1996.
10. Морева В.С. Математическое моделирование обтекания профилей с использованием новых расчетных схем метода вихревых элементов. Дис. ... канд. физ.-мат. наук. М., 2013.
11. Marchevsky I.K., Moreva V.S. Vortex element method for 2D flow simulation with tangent velocity components on airfoil surface // ECCOMAS 2012. V European Congress on Computational Methods in Applied Sciences and Engineering, e-Book Full Papers. 2012. P. 5952-5965.
12. Vaz G., Falcao de Campos J.A.C., Eca L. A numerical study on low and higher-order potential based BEM for 2D inviscid flows // Computational Mechanics. 2003. Vol. 32. Iss. 4-6. P. 327-335.
13. Макарова М.Е. Поиск аналитических решений и исследование точности расчетных схем метода вихревых элементов в двумерных стационарных задачах обтекания профилей // Инженерный журнал: наука и инновации. 2012. Вып. 4.
DOI: 10.18698/2308-6033-2012-4-166 URL: http://engjournal.ru/catalog/mathmodel/hidden/ 166.html
14. Abramowitz M., Stegun I.A. Handbook of mathematical functions with formulas, graphs, and mathematical tables. New York: Dover, 1965.
15. Лаврентьев М.А., Шабат Б.В. Методы теории функций комплексного переменного. М.: Наука, 1965. 716 с.
Кузьмина Ксения Сергеевна — аспирант, ассистент кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5).
Марчевский Илья Константинович — канд. физ.-мат. наук, доцент кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5).
Просьба ссылаться на эту статью следующим образом:
Кузьмина К.С., Марчевский И.К. Численная схема высокого порядка точности для определения интенсивности вихревого слоя при решении двумерных задач аэрогидродинамики вихревыми методами // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2016. № 6. C. 93-109. DOI: 10.18698/1812-3368-2016-6-93-109
HIGH-ORDER NUMERICAL SCHEME FOR VORTEX LAYER INTENSITY COMPUTATION IN TWO-DIMENSIONAL AEROHYDRODYNAMICS PROBLEMS SOLVED BY VORTEX ELEMENT METHOD
K.S. Kuzmina kuz-ksen-serg@yandex.ru
I.K. Marchevsky iliamarchevsky@mail.ru
Bauman Moscow State Technical University, Moscow, Russian Federation
Abstract
The study deals with the numerical simulation of two-dimensional viscous incompressible flow around airfoils by using vortex element method. The numerical scheme and the corresponding algorithm for this method usually presuppose the replacement of the airfoil with the polygon which consists of panels, and the unknown vortex layer intensity is assumed to be piecewise-constant on the panels. The accuracy of this scheme varies from O(h2) to O(h3) for different airfoils (h is the panels' length). In the present research we developed a new high-order numerical scheme. The new approach assumes the vortex layer intensity to be not piecewise-constant, but piecewise-linear or piecewise-quadratic on each panel. It is also important that the solution is not assumed to be continuous along the airfoil; it is most needed for correct simulation of the flow around airfoil with the angular points and sharp edges. Moreover, we take into account the fact that the airfoil's boundary is curvilinear: each part of the airfoil's boundary is approximated by a cubic spline instead of a straight panel. In order to obtain linear algebraic equations, we used the least squares method instead of collocation-type conditions in separate control
Keywords
Vortex method, vortex layer, high order scheme, least squares method, curvature, piecewise-polynomial approximation
points or on average on the panels. We applied higher order of accuracy Gaussian quadrature formulas for approximate integrals calculation. The results of the research show that the developed scheme has higher accuracy order than the previously known schemes. For some particular model problems (flow around circular, elliptical and Zhukovsky airfoils) this approach allows us to obtain solution with accuracy O(h5)
REFERENCES
[1] Lifanov I.K., Belotserkovsky S.M. Methods of discrete vortices. CRC, 1992.
[2] Cottet G.-H., Koumoutsakos P. Vortex methods: theory and practice. Cambridge University Press, 2000.
[3] Dynnikova G.Ya. Vortex motion in two-dimensional viscous fluid flows. Fluid Dynamics, 2003, vol. 38, no. 5, pp. 670-678.
[4] Kalugin V.T., Mordvintcev G.G., Popov V.M. Modelirovanie protsessov obtekaniya i up-ravleniya aerodinamicheskimi kharakteristikami letatel'nykh apparatov [Modeling of the flow and the aerodynamic characteristics control for aircrafts]. Moscow, MGTU im. N.E. Baumana Publ., 2011. 527 p.
[5] Andronov P.R., Guvernyuk S.V., Dynnikova G.Ya. Vikhrevye metody rascheta nes-tatsionarnykh gidrodimamicheskikh nagruzok [Vortex methods for unsteady hydrodynamic loads]. Moscow, MGU im. M.V. Lomonosova Publ., 2006.
[6] Ogami Y., Akamatsu T. Viscous flow simulation using the discrete vortex model — the diffusion velocity method. Comput. and Fluids, 1991, vol. 19, no. 3/4, pp. 433-441.
[7] Kuzmina K.S., Marchevsky I.K. On numerical schemes in 2D vortex element method for flow simulation around moving and deformable airfoils. Proc. of Summer School-Conference "Advanced Problems in Mechanics 2014". St. Petersburg, 2014. pp. 335-344.
Available at: http://www.ipme.ru/ipme/conf/APM2014/2014-PDF/2014-335.pdf
[8] Kuzmina K.S., Marchevsky I.K., Moreva V.S. On the accuracy of numerical schemes for flow simulation around airfoils with angle point. Nauka i obrazovanie. MGTU im. N.E. Baumana [Science & Education of the Bauman MSTU. Electronic Journal], 2015, no. 2, pp. 234-249. DOI: 10.7463/0215.0756954
Available at: http://technomag.neicon.ru/en/doc/756954.html
[9] Kempka S.N., Glass M.W., Peery J.S., Strickland J.H. Accuracy considerations for implementing velocity boundary conditions in vorticity formulations. SANDIA REPORT SAND96-0583 UC-700, 1996.
[10] Moreva V.S. Matematicheskoe modelirovanie obtekaniya profiley s ispol'zovaniem no-vykh raschetnykh skhem metoda vikhrevykh elementov. Diss. kand. fiz.-mat. nauk [Mathematical modeling of flow of profiles using the new calculation schemes of vortex element method. Cand. phys.-math. sci. diss.]. Moscow, 2013.
[11] Marchevsky I.K., Moreva V.S. Vortex element method for 2D flow simulation with tangent velocity components on airfoil surface. ECCOMAS 2012. V European Congress on Computational Methods in Applied Sciences and Engineering, e-Book Full Papers, 2012, pp. 5952-5965.
[12] Vaz G., Falcao de Campos J.A.C., Eca L. A numerical study on low and higher-order potential based BEM for 2D inviscid flows. Computational Mechanics, 2003, vol. 32, iss. 4-6, pp. 327-335.
[13] Makarova M.E. Calculation of flow past airfoils of simplest topology using the modified vortex-element method. Jelektr. nauchno-tekh. izd. "Inzhenernyy zhurnal: nauka i innovacii" [El. Sc.-Tech. Publ. "Eng. J.: Science and Innovation"], 2012, iss. 4.
DOI: 10.18698/2308-6033-2012-4-166
Available at: http://engjournal.ru/eng/catalog/mathmodel/hidden/166.html
[14] Abramowitz M., Stegun I.A. Handbook of mathematical functions with formulas, graphs, and mathematical tables. N.Y., Dover, 1965.
[15] Lavrent'ev M.A., Shabat B.V. Metody teorii funktsiy kompleksnogo peremennogo [Methods of the theory of functions of a complex variable]. Moscow, Nauka Publ., 1965. 716 p.
Kuzmina K.S. — post-graduate student, Assist. Professor of Applied Mathematics Department, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, Moscow, 105005 Russian Federation).
Marchevsky I.K. — Cand. Sci. (Phys.-Math.), Assoc. Professor of Applied Mathematics Department, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, Moscow, 105005 Russian Federation).
Please cite this article in English as:
Kuzmina K.S., Marchevsky I.K. High-Order Numerical Scheme for Vortex Layer Intensity Computation in Two-Dimensional Aerohydrodynamics Problems Solved by Vortex Element Method. Vestn. Mosk. Gos. Tekh. Univ. im. N.E. Baumana, Estestv. Nauki [Herald of the Bauman Moscow State Tech. Univ., Nat. Sci.], 2016, no. 6, pp. 93-109. DOI: 10.18698/1812-3368-2016-6-93-109
