Обтекание колеблющегося крыла потоком идеальной несжимаемой жидкости Текст научной статьи по специальности «Физика»

УДК 532.5

Обтекание колеблющегося крыла потоком идеальной несжимаемой жидкости

© Д.А. Крылов1, Н.И. Сидняев1, А. А. Федотов1

1 МГТУ им. Н.Э. Баумана, Москва 105005, Россия

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

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

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

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

Постановка задачи. Введем в рассмотрение неподвижную прямоугольную декартову систему координат Ох1х2 х3. Поток на беско-

нечно сти будем считать однородным с постоянным вектором скорости V«, параллельным оси х1 и направленным в положительную сторону этой оси. Пусть закон движения крыла имеет вид

хк = хк (а1,а2,1), 0 < а1 < 1, -1 < а2 < 1, (1)

где х^ и а18 (к = 1, 2,3 ; I = 1,2) — соответственно декартовы и ла-гранжевы координаты точек крыла на поверхности; I — время.

Полагаем, что бесконечно тонкое крыло схематизирует собой реальное крыло, имеющее обтекающуюся без отрыва потока закругленную кромку Ь8 и острую кромку Ь^,, с которой в поток жидкости плавно стекает вихревой след, возникающий за крылом при его движении. Поэтому считаем, что крыло является несущей поверхностью с передней кромкой (а\ = 0, -1 < а2 < 1), обтекаемой без отрыва, на которую действует подсасывающая сила [1]. С задней (а^ = 1, -1 < а2 < 1) и боковой (0 < а\ < 1, а2 = ±1) кромок в поток стекает свободная вихревая поверхность. Передняя кромка является кромкой натекания Ь8, а задняя и боковая — кромками стекания Ьп (рис. 1).

Скорости Va = vleк, индуцируемые несущей (а = s) и свободной вихревой (а = w) поверхностями в точке наблюдения с радиус-вектором r в момент времени t, вычисляют по формуле [2]

V = Va(r, t) = -L f f l Rl dal dal (2)

Здесь уа — поверхностный вектор вихря, уа = е12аПа; е12а = = |э1а хэ2а|; эа =дга/ да1а; I = 1,2 — базисные векторы лагранжевой системы координат на поверхности а; Па — вектор интенсивности завихренности поверхности а, Па = па х (У_ - У+); па — единичный вектор нормали к верхней стороне поверхности, па = (э1а х э2а) /е12а; У+, У_ — предельные значения скорости жидкости V при подходе к поверхности а с нижней («+») и верхней («-») сторон соответственно; Яа = г - га, Яа = ; г = хкек; га —

радиус-вектор точки поверхности а; га = х'аек; ек — единичные базисные векторы неподвижной декартовой системы координат, к = 1, 2, 3; аа, аа — лагранжевы координаты на поверхности а .

Пределы интегрирования для несущей и свободной вихревой поверхностей соответственно следующие: р5 = 0, ц = 1, 4 = 1,

Р* = Ьо(г), = 0, ¡* = ё, где Ь0(г) = сг; с, ё — константы, выбираемые при параметризации поверхности к исходя из удобства проведения вычислений.

Поверхностный вектор уа (а = s, *) можно разложить по базисным векторам лагранжевой системы координат:

Уа =Г1а(аа, аа,г)эа,

причем

у\ = дк/да2; у2 = -дк/да1,

где к = к(а\, а1, г) — вихревая функция, представляющая собой скачок потенциала скорости жидкости на несущей поверхности.

Особенности вихревой функции выделены в явном виде в работе [3]. Ее регулярная часть описывается гладкой функцией у/(а\, аг), определенной на всей несущей поверхности, включая ее границы.

Система уравнений для определения функций у/(а\, а%, г),

у1 а, г), а также функций х~к = х"к (а*, aw, г), задающих положение свободной вихревой поверхности, имеет вид [4, 5]

( + ук + V* - дх^/дг),к = 0, г = г^; (3)

дхкт/д! = V* + ук + Vк, г = г*0, к = 1,2,3; (4)

ду* /дг = 0, I = 1,2. (5)

Здесь п^,к = п^ек, п= п5(г^,,t); к а = 5, у — радиус-вектор точки наблюдения в случае, когда эта точка принадлежит поверхности а, к = х£о ек = хк (а^, а2о, t)ек; V» = у»ек = еь

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

хк = ху, П5 = П у, ^ 5 = ^У ; (6)

дк{а18, , t)) = -еуу(а1у, а^, t), (7)

где а1Н, а2Н, — лагранжевы координаты точки поверхности у, которая в момент времени t совпадает с точкой кромки стекания поверхности 5, имеющей лагранжевы координаты (а^, а^).

В начальный момент времени t = о необходимо задать функции у, уу и положение свободной вихревой поверхности:

у(а\, а2, о) = уо(а5, а2); хк(ау, а2, о) = хНо(аН, а2), к = 1,2,3; (8)

уУ(а2, а2,о) = у2о(а2, а2), I = 1, 2.

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

Для удобства изложения примем для лагранжевых координат поверхностей 5 и у обозначения а1, а2 и Ъ1, Ъ2 так, что а\ = а1, а2 = а2 и ау = Ъ1, ауу = Ъ2.

Правую половину несущей поверхности (о < а1 < 1, о < а2 < 1) разобъем на т х п элементов

а1 < а1 < а1+1, а2 < а2 < а1+1 (9)

сеткой координатных линий

а} = (I - 1)А1, I = 1, ..., т +1, Л1 = 1/т; а2 = С/ - 1)Л2, 7 = 1, . ., п +1, А2 = 1/п.

Сетка координатных линий свободной вихревой поверхности в текущий расчетный момент времени гц, соответствующая введенной на несущей поверхности сетке и разбивающая правую половину свободной вихревой поверхности (Ь0(гд) < Ь1 < 0, 0 < Ь2 < ё) на элементы (рис. 2)

Ь1 +1 < ь1 < Ь1, Ь2 < ь 2 < ь2+1,

определяется следующими формулами:

Ъ]г = ЪК^) = сгц, с = -1/Аг, ^ = 0,1,..., ц, ц = 1,2,...; ¿2 = Ыа2), 0 < Ь2 < п, 2 = ], ] = 1, ..., п +1, 0 < а2 < 1;

Ъ2 = ^(а1), п < Ь2 < ё, ё = т + п, 2 = т + п + 2 - /, 2 = п + 2,п + 3,...,т + п +1, I = т, т -1, ..., 1; 1 > а1 > 0.

При этом

а1 а 2

^1(а1) = + т + n, ^2 (а 2 ) = —, А1 А 2

а шаг по времени Аг считается постоянной величиной.

(10)

(11)

Рис. 2. Расчетная схема

Сетка координатных линий свободной вихревой поверхности соответствует сетке координатных линий несущей поверхности в том смысле, что координатные линии b2 = const свободной вихревой поверхно-

сти являются продолжениями координатных линий несущей поверхности (a1 = const или a2 = const). Узловые точки на задней или боковой кромке несущей поверхности являются одновременно узловыми точками передней кромки свободной вихревой поверхности, что при указанном выборе функций gi (a1) и константы c лагранжевые координаты b1 и Ь2 в узлах сетки имеют целочисленные значения.

Функция у (a1, a2, t) может быть восстановлена по своим значениям в узловых точках различными способами. В данной работе функцию y/(al, a2, t) в каждый расчетный момент времени tv будем

аппроксимировать бикубическим сплайном по переменным a1 и a2 , построенным по значениям у j функции у(a1, a2, t) в узловых точках [6]:

у, i = у j^v) = уЦ, j tn\ ii = 1, ..., m +1, ji = 1, ..., n +1, ц = 1, . ., q.

На элементе (9), который будем называть сплайновой ячейкой и обозначать индексами ij, бикубический сплайн определяется так:

у1,j(a1, a2, tM) =

m+1 n+1 ( 4 ,. ,Л( 4 , . ,Л

ZZ Z\<. 11 (*Г" ZPj (a2)

V io=1

f—2\Uo-1)

V jo =1

У Л.

(12)

'1 =1 /1=1 Здесь

а1 < а1 < аг1+1, а/ < а2 < а/+1; а1 = (а1 - а1)/А1, а2 = (а2 - а2)/Л2.

При этом учитываем симметричность течения относительно плоскости х2 = о. Последнее обстоятельство для функции у (а1, а2, t) можно записать в виде краевого условия

Иа>,о,,)+У0гй = о. (13)

Коэффициенты Я1о, ¡, ^ и Р/ /,/ зависят только от разбиения несущей поверхности на элементы (9) и могут быть найдены из решения систем линейных алгебраических уравнений порядков 4т и 4п, соответственно. Эти системы получают с использованием свойств бикубического сплайна: условий непрерывности сплайн-функции и ее производных по обеим переменным до второго порядка включительно [6].

Вначале для фиксированных значений координаты а = а

2 —

Л

(/ = 1, ..., п +1) строят частичные сплайны. По переменной а1 это будет одномерный сплайн, у которого в качестве заданных значений интерполируемой функции служат значения функции щ(а1, а2, гм) в

узловых точках координатной линии (а2 = а2). На отрезке а1 < а1 < а}+1 он имеет вид

т+1 ( 4

Щ(а1, а2, Ь) = 2 2 Я>,*, 4 (а1)'

-1\0о-1)

Л

_ _ , щ 1- (14)

1 =1 V' о =1

По переменной а2 это будет также одномерный сплайн, только заданными значениями интерполируемой одномерной функции в этом случае являются функции (14). В полосе а1 < а1 < а^, 0 < а2 < 1 он определяется выражением

п+1 ( 4

Щ(а1, а2, Ь) =2 2 Ро, >, л (а2)' Л =1 V >о=1

>о-1)

Л

Щ (а1, а2, ь). (15)

После преобразований уравнение (15) принимает вид (12).

При построении частичных сплайнов (14) и (15) для замыкания систем линейных алгебраических уравнений относительно неизвестных коэффициентов Яо, г, ^ и >, > в каждом случае необходимы еще два краевых условия, в качестве которых используют следующие [6]:

Бр"(/; хо) = Бр"(/;хм) = о;

(16)

/о+ Бр'(/; хо) = о, Бр'(/; хм) = о; /о = /(хо). (17)

Здесь Бр(/; х) — интерполяционный кубический сплайн; /(х) — интерполируемая функция; [ хо, хм ] — отрезок интерполяции, М — число частичных отрезков [ х>, х>+1] на отрезке [ хо, хм ]; х> — узловые точки. (Штрихи обозначают соответствующие производные.)

Краевые условия (17), примененные при построении сплайна вдоль оси а2 , обеспечивают симметричное обтекание крыла относительно плоскости х2 = о. Учитывая условие симметрии (13), для замыкания систем линейных алгебраических уравнений при определении коэффициентов Яо,,, ^ и Р?о, > будем использовать условия

(16) и (17) соответственно.

Запишем выражение для скорости, индуцируемой элементом (1о) свободной вихревой поверхности, в виде

! ЬЦ Ь„+1 \ э — э — Л

V., ^ „(г, О = 4- ] \ \у1 ^ + у2 ^^ (18)

4-Ь1 „„ \ ) >+1 Ь„

Из выражения (18) следует, что при известном положении в пространстве элементов поверхности . для нахождения скорости, индуцируемой свободной вихревой поверхностью, достаточно определить функции у. (Ь1, Ь2, t). В узловых точках передней кромки свободной вихревой поверхности функции у. (Ь1, Ь2, t) выразим через функцию у (а1, а2, t) [5], а в остальных узловых точках свободной вихревой поверхности — определим в соответствии с условиями (5). Положение точек свободной вихревой поверхности (функции х. (Ь1, Ь2, t) и компоненты поверхностного вектора вихря у. (функции у. (Ь1, Ь2, t)) на элементе (1о) вихревой поверхности будем искать с помощью билинейной интерполяции по значениям в узловых точках этого элемента:

х(Ь1, Ь2, 0 = хц+1,„(^)2Ц+1(Ь1)2„(Ь2) + +хц+1,„+1 (tq ) 2Ц+1 (Ь1) 2„ +1 (Ь2) + +хц„ „ (tq)2^ (Ь1)2„ (Ь2) + хц„+1 (tq)2^ (Ь1)2„+1 (Ь2 ); (19)

xц,v(tq) = х(Ь^ Ь„ , tq);

7 ь7 - Ы. 7 ь7 - Ы - -

2- (Ь7 ) = _-, 2'+1 (Ь7 ) = ^-/, I =ц,у; / =1,2,

Ь/ - Ь/ , ' Ь/ , - Ь/

II +1 I +1 I

куда вместо функции х(Ь1, Ь2, t) надо подставить функции х.(Ь1, Ь2, t), у.(Ь1,Ь2,t).

Отметим, что из интерполирования функций у. (Ь1, Ь2, t) формулами вида (19) следует, что в любой расчетный момент времени при сохранении значений функций у. (Ь1, Ь2, t) в узловых точках уравнения (5) будут выполняться в любой другой точке свободной вихревой поверхности.

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

для интегрирования которых воспользуемся методом Адамса второго порядка точности [7].

Для нахождения скорости жидкости в точке наблюдения с радиус-вектором г необходимо определить двумерные интегралы (2) по поверхностям ^ и к, которые будем вычислять приближенно, основываясь на методе ячеек [7].

Расчетный контур Ьу определим как любой замкнутый жидкий контур, начинающийся и оканчивающийся в узловой точке (о, ЪУ) задней кромки свободной вихревой поверхности (см. рис. 2). Построим конечно-разностную схему для интегрирования уравнения (7). Напомним, что если внешние массовые силы потенциальны (в нашем случае они отсутствуют), то в идеальной несжимаемой жидкости циркуляция скорости Гь по любому замкнутому жидкому контуру Ь во все время движения остается неизменной [8].

В расчетной схеме этот интегральный закон сохранения реализуем в виде условия сохранения циркуляции ГЬу по любому расчетному контуру Ьу во все моменты времени. Вычислим циркуляцию скорости по контуру Ьу в текущий tq и последующий tq+1 = tq + At расчетные моменты времени (см. рис. 2). Для этого воспользуемся теоремой Стокса. В момент времени tq в качестве поверхности, ограниченной контуром Ьу, выберем поверхность, пересекающую переднюю кромку несущей поверхности в произвольной точке (о, а 2 ), а переднюю кромку свободной вихревой поверхности — в узловой точке (Ъ\, ЬУ). Тогда

(а*,а* )

Г ьу (tq ) = | -^(а1, а2, tq )ёа1 + ^(а1, а2, tq )ёа2 + (о,а2) (о,ъУ)

+ \ -у1(Ъ\ Ъ2, ц)ёЪ1 +/1 (Ъ1, Ъ2, ц)ёЪ2, (2о)

(Ъ|,Ъ2)

где (а!, а*2) — лагранжевы координаты узловой точки на задней или боковой кромке несущей поверхности; (Ъ*, Ъ2) — то же на передней кромке свободной вихревой поверхности, совпадающей в момент времени tq с узловой точкой несущей поверхности (а^, а*2), Ъ* = Ъ1,

Ъ*2 = ЪУ = gl (а *) .

Связь между лагранжевыми координатами узловых точек (Ъ* = ЪЦ, Ъ*2 = Ъуг) и (а1 = а1, а*2 = а2) определяется формулами (11). Обозначив

(о,ЪУ)

Р = \ -у1(Ъ\ Ъ2, ^)йЪ1 +г1 (Ъ1, Ъ2, ц)ёЪ2 (Ы,ъ*2)

и учитывая, что

(О-,а* )

| -Х-2(а1, а2, tq)ёа1 + у\(а1, а2, tq)ёа2 = ^(а^, а*2, tq),

(о, а2)

перепишем выражение (2о) в виде

Гьу^ ) = к(а*, а*2, tq ) + Р. (21)

В расчетный момент времени ^+1 в качестве поверхности, ограниченной контуром Ьу, выберем поверхность, пересекающую переднюю кромку несущей поверхности в произвольной точке (о, а2), а полоску свободной вихревой поверхности, сошедшую за рассматриваемый шаг по времени, — вдоль координатной линии Ъ2 = Ъут. Тогда

ь\

Г

Lv(tq+1) = K(d,a*,tq+1)- J y2(b\b*,tq+1)db1 + F. (22)

Поскольку циркуляция по расчетному контуру Ьу одинакова во все расчетные моменты времени, то из уравнений (21) и (22) вытекает закон сохранения циркуляции по расчетному контуру Ьу в виде

ь\

к(а\, a*, tq+1) -^a1, a*, tq) = J уУ(Ь\ b*2, tq+Ц. (23)

bq+1

После подстановки первой формулы (19) (заменив предварительно функцию x(b1, b2, t) на y2(b\ b2, t)) в выражение (23) интеграл можно легко вычислить. В результате получаем уравнение

x(a*, a*2, tq+1) -к(a1, a*2, tq) =

cAtr

2

{уЗД, b2, tq) + y2(bJ+1, b2,tq+1)]^ (24)

которое представляет собой неявную конечно-разностную схему второго порядка точности [7] для уравнения (7).

Таким образом, при численном решении системы уравнений (3)-(8) с применением билинейной интерполяции вида (19) для вы-

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

Пусть в текущий момент времени tq заданы функции у(а1, а2, t),

у.(Ь1, Ь2, {) и х.(Ь1, Ь2, т. е., согласно (4) и (19), известны координаты х. (Ь1М, Ь„, tq) узловых точек свободной вихревой поверхности, значения у. (Ь1М, Ь„, tq) в этих точках и значения уь/1.

Опишем переход от расчетного момента времени tq к ^+1. Координаты узловых точек свободной вихревой поверхности получим из уравнений (4), для интегрирования которых воспользуемся методом Адамса второго порядка точности. Значения функций у.(Ь1, Ь2, t) в узловых точках, не лежащих на передней кромке свободной вихревой поверхности, определим из условий (5). Значения у.(Ь°, Ь„, tq) в узловых точках передней кромки свободной вихревой поверхности в соответствии с условиями (6) выразим через функцию у (а1, а2, t), а следовательно, через значения у /1. Последние определим из решения системы линейных алгебраических уравнений, т х п уравнений которой получим из требования выполнения условия непроницаемости несущей поверхности (3) в контрольных точках с координатами

1 (а1 + а1+1) 2 (а/2 +а/2+1) . , . ,

а0г = ' г + , а0/ = . -, г = 1, ...,m, / = 1, ...,П

2 2

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

Перепишем уравнение (24) в виде

к(а\, а*, tq+1) + ^А-уЭД+1, Ь„, tq+1) = cАt

= к(а\, а*, tq) —— у.(Ь°, Ь„, tq). (25)

Правая часть уравнения (25) известна из решения задачи в предыдущие моменты времени, а левая выражается линейно через значения у j (tq+1) . В результате получаем систему линейных алгебраических

уравнений, решение которой позволяет определить значения У+1). Поле скоростей становится известным и можно проводить дальнейшие расчеты.

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

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

Уравнения движения крыла. Движение крыла в безразмерной форме определяются соотношениями

оси вертикальных колебаний (расстояние в корневом сечении = о между точкой передней кромки крыла и точкой оси угловых колебаний положительное, если ось расположена позади передней кромки) о < Ъо < 1,5; ) = у(1) - а(1); у(1.) = аг^ (тк ът(Ш)); а(1) = = ао8т(^); ао = Ю°; о, 2 <т< 1,5; Л —удлинение крыла, Л = 4; к — амплитуда вертикальных колебаний оси угловых колебаний, о, 5 < к < 1,5; а = а(Х) — угол атаки, ао — амплитуда изменения угла атаки.

Формы крыльев в плане в плоскости х1, х2 изображены на рис. 3, а. Передняя кромка крыла 1 определяется уравнением х-. = (х|)2/3.

Опишем движение крыла в системе координат Оу1 у2у3 (рис. 3, б), в которой жидкость на бесконечности покоится и оси которой параллельны соответственно осям Ох1, Ох2, Ох3. Ось угловых колебаний, параллельная оси Оу2, движется вдоль вертикальной оси согласно закону у3 = к со8(т), а вдоль горизонтальной оси — по закону у1 = -.

х1 = а1 + Со (1 - а1 )(а2)2 - Ъо собО) + Ъо;

х2 = Л(1 - со/3)а2/2; 2

= а1 + со (1 - а1 )(а2 ) - Ъо Бт ) + к соБ(т)

(26)

Здесь о < а\ < 1; со — параметр, задающий форму крыла в плане и равный 3/4 для крыла 1 и о — для крыла 2; -1 < < 1; Ъо — положение

Крыло проектируется на плоскость у1 у3 корневой хордой АВ, а ось угловых колебаний — точкой С. Таким образом, положение крыла в каждый момент времени определяется положением оси угловых колебаний (т. е. точки С) и углом ( = ((I) наклона хорды АВ к оси Оу1. Угол между вектором мгновенной скорости V* оси угловых колебаний и осью Оу1, обозначен у = у^), а угол атаки а определяется как угол между вектором V* и хордой АВ крыла. Угол атаки а является положительным, если вектор скорости V* и нормаль п к поверхности направлены в разные стороны крыла, и отрицательным — в противном случае.

Рис. 3. Формы крыльев в плане (а) и законы движения (б): 1 — передняя кромка крыла 1; 2,3 — задняя и боковые кромки крыла 2

При переходе к безразмерным переменным за характерную длину принята длина Ь корневой хорды, за характерную скорость — скорость Уот, набегающего потока на бесконечности, за характерное время — отношение Ь / Ух. Значения параметров крыла и диапазоны их изменения установлены в результате обработки экспериментальных данных работ [10-13]. Из работы [12], в которой наиболее подробно описана кинематика хвостового плавника большого дельфина, или афалины, остается неясным положение Ьо оси угловых колебаний крыла. Вследствие этого полагали, что ее положение не зависит от времени, задается произвольно и изменяется от 0 до 1,5.

Результаты расчетов. В результате численного решения задачи определены мгновенный Ст\ и средний Ст коэффициенты силы тяги, а также и гидродинамический коэффициент полезного действия ц крыла. Коэффициент Ст^ вычисляли как отношение мгновенной силы тяги к величине рУ^Я/2, где р — плотность жидкости; З — площадь крыла. Усредненный по периоду коэффициент Ст, равен

А

у

а

б

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

Было исследовано влияние параметров к, о и Ьо на гидродинамические характеристики колеблющихся крыльев. На рис. 4, а приведены зависимости коэффициенов Ст и ] от частоты колебаний для случая, когда ось угловых колебаний проходит через точку передней кромки корневого сечения крыла (Ьо = 0), а амплитуда колебаний к = 1. При таком расположении оси как для крыла 1, так и для крыла 2 коэффициент Ст возрастает, а коэффициент ] убывает с ростом частоты колебаний. При этом коэффициент силы тяги невелик (Ст « 0,25).

Рис. 4. Зависимости коэффициентов Ст и г] от частоты колебаний о (а) и положения Ь0 оси угловых колебаний (б) для крыльев 1 (1) и 2 (2)

Существенную роль положения оси угловых колебаний для получения больших значений гидродинамических характеристик демонстрируют приведенные на рис. 4, б зависимости коэффициентов Ст и ] от параметра Ьо при к = 1 и о = 1. Коэффициент силы тяги Ст монотонно возрастает с ростом Ь0, а гидродинамический коэффициент полезного действия ] при Ьо < Ьотах монотонно возрастает, достигает максимального значения при Ьо = Ьотах и монотонно убывает при Ьо > Ьотах (Ьотах « о,9 для крыла 1 и « о, 8 — для крыла 2).

На рис. 5 приведены зависимости гидродинамических характеристик Ст и ] от частоты колебаний крыльев 1 и 2 в случае, когда ось угловых колебаний находится на задней кромке (Ьо = 1).

0,50

0 П

0,75

2,0

1,0

0,50

О 0,5 1,0 со

О 0,5 1,0 со б

а

Рис. 5. Зависимости коэффициентов Ст и ] от частоты колебаний а для крыльев 1 (а) и 2 (б) при Ь0 = 1 и И = 0,5 (1), 1,0 (2), 1,5 (5)

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

Гидродинамический коэффициент полезного действия крыла 2 (рис. 5, б) при к = 0,5 и к = 1,0 является монотонно возрастающей функцией частоты колебаний, а при к = 1,5 эта зависимость оказывается немонотонной.

На рис. 5 все зависимости Ст(а) представляют собой возрастающие функции, которые увеличиваются с ростом а тем быстрее, чем больше к. Наибольшие значения коэффициента силы тяги крыльев 1 и 2 достигаются при наибольших значениях амплитуды к и частоты а . Из анализа этих кривых следует, что колебания с большими амплитудами (к = 1,0 и 1,5) предпочтительнее колебаний с малыми амплитудами (к = 0,5) для получения больших значений коэффициента силы тяги.

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

При сравнении результатов расчетов с данными экспериментальных исследований [1о, 11] предположим, что плавание дельфинов происходит с наибольшим значением гидродинамического коэффициента полезного действия. В экспериментах [1о, 11] определяли амплитуду и частоту точки у развилки хвостового плавника дельфина. Если предположить, что движение хвостового плавника описывается уравнениями (24) с осью угловых колебаний, совпадающей с задней кромкой, то эксперименты [1о, 11] дают амплитуду и частоту оси угловых колебаний моделей хвостового плавника. Предположение о положении оси угловых колебаний является естественным, поскольку из приведенных на рис. 4 и 5 расчетных данных следует, что для плавания с наибольшим значением коэффициента ] ось угловых колебаний должна располагаться вблизи задней кромки. В указанных экспериментах в основном наблюдали амплитуду и частоту колебаний точки у развилки хвостового плавника, соответствующие к » 1, о; со» 1,1 [1о] и о» о, 9 [11].

Расчеты, проиллюстрированные на рис. 5, а, показывают, что в случае Ьо = 1 наибольшее значение коэффициента ] достигается при к = 1,о и о » 1,15. Таким образом, результаты расчетов для крыла 1 находятся в согласии с предположением о том, что дельфины при длительном равномерном движении (режим равномерного прямолинейного движения [4]) используют наиболее рациональный механизм плавания. В работах [1о, 11] указано, что с увеличением амплитуды колебаний хвостового плавника у дельфинов наблюдается тенденция к уменьшению частоты колебаний. Согласно графикам на рис. 5, а, с позиции гидродинамики это происходит по той причине, что при увеличении амплитуды колебаний наибольший коэффициент ] достигается при более низких частотах колебаний.

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

Большой интерес представляет исследование формирования и движения свободной вихревой поверхности за колеблющимся кры-

лом. Поскольку такое исследование усложняется трехмерностью картины обтекания, ограничимся рассмотрением свободной вихревой поверхности в плоскости симметрии течения [14, 15]. На рис. 6 показано положение свободной вихревой поверхности в плоскости х1 х3 при г = 3^ за крылом 1, закон движения которого определяется уравнениями (24) со следующими значениями параметров: к = 1, Ьо = 1, а = 1. Здесь отрезком АВ изображена корневая хорда крыла (точка А соответствует кромке натекания корневой хорды), точками показаны положения узловых точек (рядом с точками приведены их ла-гранжевы координаты Ь/ = -Ь/ , / = 0, 1, 2, ..., д). Анализ результатов

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

в момент времени т = 3я

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

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

Работа выполнена при поддержке РФФИ (проект № 13-08-00538-а).

ЛИТЕРАТУРА

1. Краснов Н.Ф. Аэродинамика. Ч. I. Москва, Высш. шк., 1980, 495 с.

2. Зайцев А.А. Теория несущей поверхности. Математическая модель, численный метод, расчет машущего полета. Москва, Наука, Физматлит, 1995, 160 с.

3. Chopra M.G., Kambe T. Hydromechanics of Lunate-Tail Swimming Propulsion. Part 2. J. FluidMech, 1977, vol. 79, no. 1, pp. 49-69.

4. Федотов А.А. Эффективность работы хвостового плавникового движителя. Докл. АН СССР, 1987, т. 293, № 1, с. 48-51.

5. Федотов А.А. Структура вихревого следа за крылом, работающим в режиме нормального трепещущего полета. Вестник МГУ. Серия 1. Математика. Механика. 1990. № 3. С.42-46.

6. Завьялов Ю.С., Квасов Б.И., Мирошниченко В.А. Методы сплайн-функций. Москва, Наука, 1980, 352 с.

7. Калиткин Н.Н. Численные методы. Москва, Наука, 1978, 512 с.

8. Седов Л.И. Механика сплошной среды. Т. I. Москва, Наука, 1976, 536 с.

9. Белоцерковский С.М., Ништ М.И. Отрывное и безотрывное обтекание тонких крыльев идеальной жидкостью. Москва, Наука, 1978, 351 с.

10. Пятецкий В.Е., Каян В.П. О кинематике плавания дельфина-афалины. Бионика, 1975, вып. 9, с. 41-46.

11. Каян В.П., Пятецкий В.Е. Кинематика плавания дельфина-афалины в зависимости от режима ускорения. Бионика, 1977, вып. 11, с. 36-41.

12. Каян В.П. О гидродинамических характеристиках плавникового движителя дельфина. Бионика, 1979, вып. 13, с. 9-15.

13. Козлов Л.Ф. Теоретическая биогидродинамика. Киев, Вища шк., 1983, 238 с.

14. Федотов А.А. Расчет вихревой структуры за крылом, работающим в режиме создания силы тяги. Альманах современной науки и образования. Тамбов, 2008, № 7(14), с. 225-229.

15. Федотов А.А. Исследование вихревой структуры за колеблющимся крылом. Современные методы теории краевых задач, материалы Воронежской весенней математической школы «Понтрягинские чтения - XIX». Воронеж, ВГУ, 2008, с. 213-214.

Статья поступила в редакцию 21.02.2013

Ссылку на эту статью просим оформлять следующим образом:

Д.А. Крылов, Н.И. Сидняев, А.А. Федотов. Обтекание колеблющегося крыла потоком идеальной несжимаемой жидкости. Инженерный журнал: наука и инновации, 2013, вып. 2.

URL: http://engjoumal .ru/catalog/appmath/hidden/607.html

Крылов Дмитрий Александрович — ассистент кафедры «Высшая математика» МГТУ им. Н.Э. Баумана; аспирант дневного отделения; автор 14 научных работ; сфера научных интересов: численные методы, математическое моделирование. e-mail: dmitrykrylov@rambler.ru

Сидняев Николай Иванович — д. т. н., проф., заведующий кафедрой «Высшая математика» МГТУ им. Н.Э. Баумана; автор 200 научных работ; сфера научных интересов: численные методы, уравнения математической физики, механика жидкости, газа и плазмы, вероятность и статистика. e-mail: sidnyaev@yandex.ru

Федотов Анатолий Александрович — доцент кафедры «Высшая математика» МГТУ им. Н.Э. Баумана; канд. физ.-мат. наук; автор более 50 научных работ; сфера научных интересов: гидродинамика крыла, численные методы, уравнения математической физики. e-mail: le-tail@list.ru