CC BY
89
Том XXXIII
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
_____
М3—4
УДК 533.6.013.2.011.32:629.7.025.73 629.735.33.015.3.025.73
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВЛИЯНИЯ СХОДА ВИХРЕЙ НА НЕСТАЦИОНАРНЫЕ АЭРОДИНАМИЧЕСКИЕ ХАРАКТЕРИСТИКИ ПРОФИЛЯ ПРИ ЕГО ПРОИЗВОЛЬНОМ ДВИЖЕНИИ
А. Н. ХРАБРОВ
Рассматривается произвольное неустановившееся движение тонкого профиля на малых углах атаки. С задней кромки профиля сходит вихревой след, поэтому нестационарные аэродинамические нагрузки зависят не только от мгновенных значений кинематических параметров движения, но и от всей предыстории движения. Для вычисления нестационарных аэродинамических характеристик разработано математическое обеспечение, основанное на решении интегрального сингулярного уравнения Вольтерра первого рода для нахождения функции разрыва скоростей в следе за профилем.
На основании полученных точных численных решений показано, что в задачах динамики полета и аэроупругости влияние вихревого следа на аэродинамические нагрузки может моделироваться решением системы обыкновенных дифференциальных уравнений для дополнительных внутренних переменных. Точность моделирования определяется порядком дополнительной динамической системы. Представлены результаты моделирования в Сравнении -г с точным численным решением как для квазиступенчатого изменения кинематических параметров движения, так и для более сложного квазислучайного движения с широким спектром.
Классическая теория нестационарного движения профиля в идеальной несжимаемой жидкости была сформулирована в работах Бирнбаума, Кюсснера (немецкая школа), Лаврентьева, Келдыша, Седова (московская школа), Сирса, Кармана, Теодорсена (американская школа). Подробный обзор работ и ссылки приведены в монографии [1]. Были решены задачи о нахождении нестационарных аэродинамических нагрузок при малых гармонических колебаниях профиля с различными частотами, а также при мгновенном изменении угла атаки. Показано, что вихревой след, сходящий с задней кромки профиля при неустановившемся движении, существенным образом влияет на аэродинамические нагрузки и хранит в себе всю предысторию движения. В настоящей статье представлено разработанное математическое обеспечение, позволяющее решать сингулярное интегральное уравнение Вольтерра первого рода [2] для функции разрыва скорости в следе и находить аэродинамические нагрузки при произвольном движении профиля. В основе метода лежит математический алгоритм с использованием быстрого преобразования Фурье [3]. В статье кратко описана постановка задачи и приведены примеры расчетов.
В задачах динамики полета и аэроупругости необходимо иметь простые математические модели нестационарной аэродинамики, которые можно было бы применять на каждом шаге интегрирования по времени основной динамической системы. Широко используемые в настоящее время модели с вращательными и нестационарными аэродинамическими производными не позволяют учесть влияние предыстории движения. Использование переходных функций или решение интегрального уравнения на каждом шаге решения динамической задачи затруднительно ввиду серьезного увеличения времени расчета. В работах [4], [5] предлагается использовать для описания изменения нестационарных аэродинамических нагрузок дополнительные внутренние переменные. В настоящей работе проведено подробное исследование такого подхода, при этом
результаты моделирования с разным порядком аппроксимации сравниваются с полученными точными численными решениями интегрального уравнения.
1. Рассмотрим неустановившееся движение тонкого симметричного профиля в несжимаемой жидкости без отрыва потока. Схема течения представлена на рис. 1. Используем две системы координат — подвижную Оху, которая связана с профилем, и неподвижную Ох\у\. Подвижная система координат движется относительно неподвижной с постоянной скоростью основного движения ?70. Для обезразмеривания задачи считается, что хорда профиля равна единице и скорость основного движения и0=1. В этом случае безразмерное время t измеряется в пройденных хордах профиля. Начало подвижной системы координат связано с условным центром тяжести профиля, находящимся на расстояшш хо от носка профиля. Вращение профиля по т&нгажу происходит также относительно этой точки. Таким образом, в подвижной системе координата передней кромки профиля х0, задней кромки х0-1. Точка Ь0 в подвижной системе координат соответствует координате, в которой задняя кромка профиля находилась в момент начала неустано-вившегося движения. Неустановившееся движение профиля как твердого тела характеризуется зависимостями от времени его угла атаки а(?) и безразмерной угловой скорости вращения по тангажу Г2(0- Считается, что аиО малы (а«1, 0«1) и профиль тонкий. Вследствие этого граничные условия на профиле можно снести на отрезок (*0-1, *о) оси Ох, а граничные условия на вихревой пелене, сходящей с профиля, — на отрезок (Ь0, х0-1).
Решение задачи о вычислении нестационарных аэродинамических нагрузок, действующих на профиль при его произвольном неустановившемся движении, при сделанных выше предположениях проводится аналогично работе [2], в которой получены выражения для нестационарных аэродинамических нагрузок для случая хо=1/2. В комплексной плоскости г=хНу находится функция скоростей абсолютного течения жидкости
— = м-ги. сЬ
Эта функция регулярна вне отрезка (Ь0, х0), исчезает на бесконечности — ->0 и удовле-
ск
творяет следующим граничным условиям. При переходе через линию разрыва скоростей (вихревую пелену) на отрезке (й0, 1-л0) вертикальная составляющая скорости V непрерывна. Скорость и конечна в точке ;со-1 (условие Чаплыгина — Жуковкого в задней кромке). На обеих сторонах отрезка (*о-1, *о) вертикальная скорость одинакова и имеет заданное значение у-ь„(х, I), определяемое движением крыла. Таким образом, задачу нахождения функции — во всей плоскости
ск
можно свести к ее определению в верхней или нижней полуплоскости с добавочным граничным условием
и = 0, ле(-оо,60)и(л:0,оо).
Введем обозначения щ, иг для значений продольной скорости при подходе к разрезу (Ь0, *о) снизу и сверху. С учетом предыдущего условия имеем
и2(х) = -щ(х).
Пусть Г(0 — циркуляция по контуру, охватывающему профиль с частью вихревой пелены и пересекающему вихревую пелену в точке с координатой 5 в неподвижной системе координат. Тогда имеем
ds
= 2u2(s) .
Для нахождения функции <^/<к воспользуемся формулой Коши с учетом необходимых особенностей на передней и задней кромках профиля:
dw _ 1 \z-xQ + 1 dz m'J z — jc0
*o
,*o-i
s0
“гО)
^-гЦ-XQ+l
■ds
(1)
Здесь 5о и ^1 — координаты начала и конца вихревой пелены в неподвижной системе координат, а переменные £ и 51 связаны соотношением:
(2)
Таким образом, имеем аналитическое выражение для функции скоростей абсолютного движения жидкости в зависимости от нормальной скорости потока на профиле (известная функция при заданных законах изменения а(г) и £1(0 и разрыва скорости в вихревой пелене за профилем u2(s)). Эта последняя функция нуждается в дополнительном определении.
Для нахождения функции разрыва скорости в следе и вычисления аэродинамических нагрузок понадобится разложение функции скоростей в ряд Лорана в окрестности бесконечно удаленной точки. Так как на бесконечности — -> 0, то искомое разложение имеет вид:
dz
dw _ Г0 1 , ic2 , ic3 , dz 2m z z1 z3
коэффициенты которого Го, C2 и Сз легко вычисляются:
(3)
с, =-
2я
Го =-2
2, I h-^+1
■ds
ds
c3 -
8я
*o
!
*0-1
f vn [4x0 -1 + 4^(1 + 2^)]J t d\+
Il-Xn + £
I
+ JM2[4x0-l + 4^(l + 2O]J-
so ■'
Xq+1
ds
(4)
Величина Г0 представляет собой циркуляцию по бесконечно удаленному контуру, которая должна быть постоянна во времени. В момент начала неустановившегося движения, когда следа еще нет, циркуляция по профилю может быть выражена в виде соотношения:
г»=“2]"-«Р^Т7^'
(5)
Здесь vnQ — значение вертикальной скорости на профиле в момент начала неустановивше-гося движения. Поэтому первое выражение из (4) может быть переписано в виде:
(6)
где в левой части стоит интеграл от неизвестной функции разрыва скоростей в вихревом следе, а в правой части — интеграл от возмущений нормальной скорости на профиле по сравнению со скоростью на профиле, существовавшей при установившемся движении. При заданном законе неустановившегося движения профиля правая часть выражения (6) — известная функция. Таким образом, соотношение (6) является интегральным уравнением для нахождения неизвестной функции и2(я). Если начало неподвижной системы координат выбрать в точке, совпадающей с положением задней кромки профиля в момент начала неустановившегося движения (50 = О, = г), то с учетом связи (2) уравнение (6) может быть переписано в виде:
Данное уравнение является интегральным уравнением Вольтерра первого рода со слабой (интегрируемой) сингулярностью. Ниже будет представлена программа решения такого интегрального уравнения.
Так как известно разложение функции — —>0 в окрестности бесконечно удаленной точ-
сЬ
ки (3), формулы для коэффициентов подъемной силы и момента тангажа относительно начала подвижной системы координат могут быть записаны в виде [2]:
/ + 1
(7)
Достаточно громоздкие вычисления приводят к следующим результатам для коэффициента подъемной силы:
•*0-1 1 *о—1
и коэффициента момента тангажа относительно начала подвижной системы координат
4->*
Хо-1
-2Ь0 -^| ] -^(*оЧХЬ*о+М + *0-1
2
г
\
Видно, что нестационарные аэродинамические нагрузки, действующие на профиль, имеют три составляющие:
Здесь квазистационарные аэродинамические нагрузки суХ и тл определяются мгновенным
распределением скорости потока вдоль поверхности профиля. Эта часть аэродинамических нагрузок выражается в формулах (8) и (9) интегралами от нормальной скорости и„(х, г) на профиле. Так называемые присоединенные массы суг и тл выражаются через интегралы от производной по
времени Собственно нестационарные аэродинамические нагрузки су3 и тл зависят от раз-<М
рыва скоростей в следе (интегралы по следу от функции м2)- Именно эти члены дают зависимость аэродинамических нагрузок от всей предысторий движения.
В случае заданной формы профиля и известного закона движения первые две составляющие аэродинамических нагрузок могут быть вычислены явно. Для вычисления третьей нестационарной составляющей необходимо сначала решить интегральное уравнение Вольтерра (7) с заданной функцией /(1).
Произвольное движение профиля описывается законами изменения по времени угла атаки а(Г) и безразмерной угловой скорости тангажа О(Г). При этом распределение нормальной скорости для пластинки в каждый момент времени имеет следующий вид (считается, что вращение происходит относительно начала координат):
С использованием этих выражений могут быть вычислены определенные интегралы по поверхности профиля, входящие в выражения (8) и (9). В этом случае из (10) получим следующие выражения:
°У ~ + Су2 + су3 ’
(10)
т2 - тл + тг2 + тгу
р„(;с,0 = -а(0 + С2(0* •
Соответствующая производная по времени может быть выражена следующим образом:
Л
су (0 = Су а(0 + с^П(0 + Cy6.it) + с“П(г) + су3, т2(0 = /и“а(0 + трП(0 + т* а(0 + т^О(0 + [*о
(П)
Здесь введены обозначения для квазистационарных аэродинамических производных:
Аэродинамические нагрузки, вызванные присоединенными массами жидкости, описываются производными с точками:
Для членов, описывающих влияние вихревого следа, с учетом замены переменных (2) введено обозначение:
(12)
Таким образом, для полного вычисления аэродинамических нагрузок осталось решить интегральное уравнение (7) и найти функцию щ^).
2. В работе [3] приведен алгоритм численного решения скалярного интегрального уравнения Вольтерра со слабой сингулярностью
Предполагается, что ядро к(О, а такжеf(t)v^g{t,y) — непрерывные функции аргумента Г при />0. Тогда решение у(1) также непрерывная функция. Для вычисления свертки используется быстрое преобразование Фурье, поэтому задача решается на равномерной сетке с шагом А, число
Функция /(* ) вычислена для движения профиля как твердого тела, при этом считается, что неустановившееся движение начинается из стационарного при а=0, 0=0. Это предположение не накладывает дополнительных ограничений, так как можно рассматривать получающееся решение линейной задачи как возмущение соответствующего стационарного решения.
На основе алгоритма [3] была разработана программа в пакете МАТЬАВ для численного решения рассматриваемого интегрального уравнения.
Для проверки рассматриваемого алгоритма решения интегральных уравнений Вольтерра первого рода было рассмотрено интегральное уравнение Абеля
Сравнение численного решения данного уравнения, полученное на сетке ти=10 (1031 узел) для интервала 0 < £ < 8, с соответствующим точным решением показало .прекрасное совпадение результатов.
(13)
узлов которой равно 2т+7. Исследуются условия сходимости и устойчивости. Показывается, что скорость сходимости решения к точному составляет 0(й4). Там же предлагается модификация этого метода для решения аналогичной задачи первого рода:
о
при этом дополнительно предполагается, что к(0) * 0 и дg/дy * 0. Уравнение (13) совпадает с уравнением (7), если положить
к(* - 5) = У^-5 + 1,
?(^7(‘«)) = Ы2('уХ
имеющее точное решение
У(0=-Т=£-
Кроме этого было численно найдено решение уравнения
1 '
имеющее точное решение
хо=г
Особенностью последнего решения является его сингулярность в окрестности точки /=0. Численно решение находилось на интервале 0<t<2 при /я=10 (ТУ=1031) и /и =12 (N=4103). Отличие численного решения от точного заметно только в малой окрестности сингулярности. В остальном интервале результаты практически совпадают. В окрестности же сингулярности решения в численных решениях возникает некоторая колебательная неустойчивость, характерная для интегральных уравнений первого рода. Чем выше порядок сетки, тем окрестность сингулярности, где возникает колебательность решения, меньше. Уменьшается с измельчением сетки и амплитуда колебательности.
Время получения решения на компьютере Репйит 133 МГц составляет 1—2 минуты при сетке порядка /и =10. С увеличением порядка сетки (измельчением шага) время решения возрастает нелинейным образом.
Результатом решения интегрального уравнения (7) является непрерывная функция ), определяющая скачок скорости при пересечении следа. Для нахождения нестационарных составляющих аэродинамических нагрузок в соответствии с выражениями (10) и (12) необходимо вычисление определенного интеграла по вихревому следу. Этот интеграл также является слабо сингулярным, для его вычисления воспользуемся равенством
( /—А /
Г и2{з)(к Г «2 (■?)*& Г Ц2 (■?)<&
о^-я)(*-я + 1) о -■* + !)
где А — шаг сетки, на котором получено решение и2^). Первый интеграл особенностей не содержит и может быть легко вычислен стандартными численными методами. Особенность второго интеграла интегрируема, для его вычисления воспользуемся заменой переменных t-s = x, тогда имеем:
г и2(з)с1з _ [и2(1. -Х)сЬс -■*)(*-■* + !) І Й> + 1)
(14)
С точностью до 0(И2) можно записать, что и2 (1-х) = и2 (?) - хс1и2 /Ш и 1/VI + х = 1-х/2 . Подставляя эти соотношения в выражение (14) при сохранении того же порядка точности 0(И2), получим
1 Ул 2>1 1 Ал г)
Эти выражения используются в дальнейшем при вычислении интеграла по следу для расчета нестационарной составляющей суз-
Воспользуемся предлагаемой методикой для расчета нестационарных аэродинамических характеристик при произвольном законе изменения угла атаки. В качестве примера выбрана следующая функция изменения угла атаки:
Рис. 2. Функция разрыва скорости за профилем при квазиступенчатом изменении угла
атаки с различным темпом
Эта непрерывная функция представляет собой синусоидальное увеличение угла атаки от а=0 при * < /) до а=1 при В качестве максимального значения угла атаки была выбрана единица, так как задача линейна и решение для другого максимального угла атаки может быть найдено просто умножением полученного решения на соответствующий коэффициент. Данная ква-зиступенчатая функция имеет непрерывную первую производную.
При заданном законе движения профиля сначала решалось интегральное уравнение (7) и находилась функция иг(0- Затем с помощью выражений (11) находились нестационарные аэродинамические нагрузки.
На рис. 2 представлены решения интегрального уравнения для /1=0,1 и Д/=*2-/1=1, 0,5 и 0,25|, что соответствует увеличению угла атаки в течение времени, за которое профиль проходит одну свою хорду, половину хорды и четверть хорды. Решение интегрального уравнения проводилось до / = 4. Сетка решения выбиралась порядка ти=10, что соответствует числу расчетных точек 1031. Видно, что с увеличением темпа изменения угла атаки разрыв скоростей в следе (завихренность) существенно возрастает.
На рис. 3 в верхней части показаны полученные при А/ = 0,5 зависимости суммарного коэффициента подъемной СИЛЫ Су и отдельных его составляющих Су], Су2 И Суз- В нижней части этого рисунка показаны соответствующие зависимости для суммарного коэффициента момента тангажа т2 и его составляющих тг\, тл и /иг3 относительно начала координат, расположенного в точке *0=0,3. Видно существенное запаздывание развития суммарных аэродинамических нагрузок по сравнению с квазистационарным значением, обусловленное влиянием вихревого следа,
Анализ этих, и других полученных результатов для разных темпов увеличения угла атаки показывает, что при, увеличении а численные решения для коэффициента подъемной силы стремятся к известному теоретическому результату для ступенчатого изменения угла атаки (функция Вагнера [1]). Из теории следует, что в несжимаемой жидкости при мгновенном возрастании угла атаки в первоначальный момент времени подъемная сила будет равна бесконечности. В следующий момент времени коэффициент подъемной силы будет равен половине, а затем постепенно стремится к стационарному значению, соответствующему новому углу атаки.
В данных расчетах полагалось, что угловая скорость профиля цо тангажу 0 = 0. Проводились также численные расчеты и для случая квазиступенчатого изменения 0(0 при постоянном
ю
в в 4
3 2 1 О
-1
-2
■*0 0,8 1 1,5 2 V 3 *
УМ г V
ч
о
■V -V ■V
-10 1 1,5 2 2£ 3 3,8 4
Рис. 3. Составляющие аэродинамических коэффициентов подъемной силы и момента тангажа при квазиступенчатом изменении угла атаки за время Д/ = 0,5
а = 0. Результаты качественно похожи на данные, полученные при изменении угла атаки. Вихревой след также существенно влияет на суммарные значения аэродинамических нагрузок.
3. Для разработки простой математической модели, позволяющей описывать влияние вихревого следа на нестационарные аэродинамические характеристики и удобной для использования в задачах динамики, рассмотрим классическое решение нестационарной аэродинамики об установившихся колебаниях профиля с малой амплитудой при нулевом угле атаки с различными частотами. Известно [1], [2], что в этом случае нестационарная составляющая коэффициента подъемной силы выражается следующим образом:
су3 =[С(ю)-1]
2яа + 2л| —-дс0.
(15)
Здесь С(со) — комплекснозначная функция Теодорсена, зависящая от безразмерной частоты колебаний. Она выражает амплитуду и сдвиг фаз нестационарной аэродинамической составляющей относительно соответствующего квазистационарного значения. В теории получено ее выражение через специальные функции, вычислены предельные значения при очень малых (со->0) и очень больших (©—»со) частотах колебаний профиля
(1«)
При больших частотах влияние вихревого следа за профилем особенно значительно и приводит к снижению в два раза производной коэффициента подъемной силы по углу атаки и угловой скорости тангажа.
В связи со сказанным выше функцию Теодорсена можно рассматривать как некоторую передаточную функцию для нестационарной составляющей аэродинамических нагрузок.
Будем аппроксимировать функцию Теодорсена дробью с полиномами порядка п в числителе и знаменателе
\ + а1р + а2р2 + ... + апр" 1 + \р + Ъ2р2 +... + 2 апр"
ся(р) = ; у _, 2 о7)
где р — оператор Лапласа. При записи выражения (17) были учтены предельные выражения (16).
Задача состоит в определении необходимого порядка аппроксимации и и нахождения неизвестных коэффициентов а, и Ъг. Для ее решения можно сформулировать соответствующую задачу минимизации. Функционал, требующий минимизации, можно составить из суммы квадратов модулей рассогласования точной и аппроксимированной функции Теодорсена:
Ф = ]Г|С(шу)-СЛсоу)|2. (18)
Для вычисления этого функционала сумма бралась по значениям частоты щу'= 1,..., 500, равномерно распределенным в интервале (0,10). При минимизации использовался симплексный метод Нелдера — Мида прямого поиска.
Для первого порядка аппроксимации имеем
С1(р)=^^- = - + -Ы- = --^£-, (19)
1 1 + 2тр 2 1 + Тр р + р!
т. е. при минимизации необходимо определить только одну неизвестную константу Т= 2т. Результаты расчетов показали, что Т = 2,5010 (Р1 = 1/Т), при этом функционал рассогласования равен Ф = 0,4513.
При нахождении аппроксимации второго порядка выражение (17) также было разложено на элементарные слагаемые, так как нахождение неизвестных параметров для вида (17) затруднительно в связи с возникновением ряда локальных экстремумов. При приведений же к виду
(.(МЛ-------^----1 (20)
2 \ + ТуР 1 + Т2р р + Р! р + р2
расчеты вне зависимости от начальных условий приводят к минимальной точке с координатами а = 0,2211, Т\ = 1,1631 и 7г = 5,9771 ({3] = 1 /Т\, Рг = 1/Гг). При этом минимум функционала составляет Ф= 0,0184.
Также был рассмотрен случай аппроксимации функции Теодорсена полиномами третьего порядка. В этом случае количество неизвестных параметров равно пяти, и при разложении на элементарные слагаемые аппроксимация имеет вид:
ад-1+^_+_»_+4*=И±а.
3 2 \Щр 1 + Т2р 1 + Тгр
Ар Вр (0,5 -А-В)р
р + (31 р + р2 р + Рз
В этом случае минимизация функционала рассогласований приводит к решению А = 0,0936, 5 = 0,2915, Тх = 0,7331, Т2 = 2,6330, Г3 = 13,8170 (^ = 1/Гь р2 = 1/Т2, 0з=1/Г3). Значение функционала в точке минимума составляет Ф=0,0015. Видно, что с повышением порядка аппроксимации функционал рассогласования значительно уменьшается.
В верхней части рис. 4 показано сравнение результатов аппроксимации с точными значениями раздельно для действительной и мнимой частей функции Теодорсена. Видно, что аппроксимация первого порядка справедлива только для использования при качественном моделирова-
Рис. 4. Аппроксимация разного порядка для функций Теодорсена и Вагнера: разными линиями показаны результаты различных аппроксимаций; Exact — точное
решение
нии. Аппроксимации второго и третьего порядка позволяют говорить и о количественном совпадении.
Полученные аппроксимации функции Теодорсена можно трансформировать в аппроксимации функции Вагнера, второй классической функции нестационарной теории профиля. Функция Вагнера описывает, как уже упоминалось выше, реакцию подъемной силы профиля на мгновенное увеличение угла атаки.
Регулярная часть функции Вагнера является переходной функцией для коэффициента подъемной силы. Используя связь переходной функции с интегралом Фурье, можно записать
00
с“(©) = е> |щт)8шют</т, (22)
о
где с“(со) — реакция коэффициента подъемной силы при малых колебаниях угла атаки с частотой <о. Из полученных выше соотношений для коэффициента подъемной силы можно записать
Су (о) = 2 л Ле С(ю), (23)
или с учетом аппроксимации первого порядка (19) получим
г
\У(0 = 4
\
2т
Л Л
-----е
2 4
V
= 2л
1—е 2
1 -IV
что является аппроксимацией функции Вагнера первого порядка. Выполняя аналогичные вычисления для аппроксимации второго порядка (20), получим
У/Ц) = 2п
1 -ае ^ - Г'1 ^ — а е~Р2‘
12 ,1
для третьего порядка (21) соответственно
1Ы(() = 2т1
1-Ае^1 -Вё~^ - (1 } --А —В е-Рз»
В нижней части рис. 4 показаны результаты сопоставления точного выражения для функции Вагнера с ее аппроксимациями первого, второго и третьего порядка. Таким образом, рассматривая функцию Теодорсена как передаточную функцию динамической системы, представляющей собой нестационарные нагрузки на профиле в несжимаемой жидкости с учетом влияния вихревого следа, можно получить конечномерные аппроксимации, описывающие полученные результаты с необходимой точностью. Причем, эти аппроксимации в частотной области приводят к соответствующим правильным аппроксимациям и во временной области.
В соответствии с выражением (15) и аппроксимациями функции Теодорсена (19) — (21) для моделирования влияния вихревого следа за профилем будем иметь следующие системы дифференциальных уравнений:
;с + р1л: = -л<х-л су з -х
—*о
О,
для первого порядка аппроксимации;
для второго порядка аппроксимации;
(Ъ \ • х2 + $2х2 = -2Впа - 2Вп — х0 О,
И
су 3=х1+х2+х3
для третьего порядка аппроксимации.
Следовательно, описание влияния предыстории движения на нестационарные аэродинамические характеристики сводится к введению дополнительных динамических внутренних переменных. Характерные времена этих переменных определяются при аппроксимации функции Теодорсена. Эти переменные теоретически с одинаковым запаздыванием влияют на подъемную силу и момент тангажа профиля.
На рис. 5 показаны результаты точного численного расчета для коэффициента су3 на основе решения интегрального уравнения Вольтерра для квазиступенчатого изменения угла атаки с различным темпом, а также результаты моделирования с помощью систем дифференциальных уравнений (24)—(26). Видно, что при моделировании нестационарного члена системой второго и третьего порядка результаты получаются уже достаточно близкими друг к другу и к точному решению.
суЪ
°#Г
о
о
1
2
3
Рис. 5. Различные порядки моделирования нестационарного члена при квазиступенча-том изменении угла атаки (обозначения те же, что на рис. 4)
а(0, «2(0
Рис. 6. Моделирование нестационарного члена при квазислучайном изменении угла атаки (обозначения те же, что на рис. 4)
Для демонстрации того, что предлагаемая математическая модель нестационарных аэродинамических характеристик работает не только для квазиступенчатого изменения кинематических параметров движения, но и для более сложных зависимостей от времени, рассмотрим изменение угла атаки а(/) в виде ряда Фурье с ограниченным спектром
N / >1/2
а^=3т
Аг=1 v z
( 2 nkt л cosl'— + е* |>
(27)
где рк — относительная мощность сигнала для &-й гармоники рк = 1 ^, а 0* — ее фазовый
сдвиг. В работе [6] решена задача о выборе зависимости фазового Сдвига гармоники от ее номе-
ра для того, чтобы минимизировать автокорреляционную функцию сигнала. Этим достигается генерация гладкого квазислучайного сигнала. Для равномерного спектра (pk=\/N) необходимо задать
пк2 (28)
Здесь 01 — сдвиг фазы первой гармоники. В настоящей работе фаза 0j выбиралась таким образом, чтобы обеспечить для зависимости угла атаки от времени начальное условие а = 0.
Полученная в соответствии с формулами (27) и (28) зависимость a(t) показана сплошной линией на верхнем графике рис. 6. При этом было выбрано Т- 60 (профиль проходит 60 хорд) и N= Ю, т. е. сигнал имеет плоский спектр гармонических составляющих в диапазоне от минимальной безразмерной частоты 0,105 до максимальной частоты 1,047. На этом же графике пунктирной линией показано решение интегрального уравнения Вольтерра для функции разрыва скорости в следе за профилем мг(0 для рассматриваемого закона изменения угла атаки. Данное решение получено на сетке порядка m =10 с числом ячеек 1031. На нижнем графике этого же рисунка показана нестационарная составляющая коэффициента подъемной силы су3, описывающая влияние вихревого следа. Различными видами линий показаны результаты моделирования с различным порядком аппроксимации и точного численного решения задачи. Видно, что результаты, полученные с использованием модели второго порядка, практически совпадают с точным решением. Проводились также расчеты для более широкого диапазона частот квазислучайного сигнала с теми же качественными результатами.
Таким образом, показано, что влияние предыстории движения, обусловленное наличием вихревого следа за профилем, на нестационарные аэродинамические нагрузки при его произвольном движении может моделироваться с необходимой точностью решением системы обыкновенных дифференциальных уравнений. Система второго порядка с двумя дополнительными внутренними переменными демонстрирует хорошее совпадение результатов с полученными данными точного численного решения задачи на основе интегрального уравнения Вольтерра при различных видах движения как при квазиступенчатом изменении кинематических параметров, так и при квазистохастическом движении с широким спектром.
Работа выполнена при поддержке Российского фонда фундаментальных исследований (проект 99—01—00042).
ЛИТЕРАТУРА
1. Некрасов А. И. Теория крыла в нестационарном потоке. — М.—Л.: Изд. Академии наук СССР.—1947.
2. С е д о в JI. И. Плоские задачи гидродинамики и аэродинамики. Изд. третье, переработанное. — М.: Наука. — 1980.
3. Hairer Е., Lubich Ch., Schlichte М. Fast numerical solution of weakly singular Volterra integral equations // J. of Computational and Applied Mathematics. —1988. Vol. 23.
4. К a r p e 1 M. Design for active flutter suppression and gust alleviation using stale-space , ■ • aeroelastic modeling//J. of Aircraft. — 1982. Vol. 19, N 3.
5.Leishman J. G., Nguyen K. Q. State-space representation of unsteady airfoil behavior//AIAA J. —1985. Vol. 28, N 5.
6. Schroeder M. R. Synthesis of low-peak-factor signals and binary sequences with low autocorrelation // IEEE Transaction on Information Theory. — January, 1970.
Рукопись поступила 20/VIII2001 г.
