Нестационарные аэродинамические нагрузки на профиле в произвольном вертикальном порыве и их моделирование с помощью обыкновенных дифференциальных уравнений Текст научной статьи по специальности «Математика»

Том XXXVIII

УЧЕНЫЕ ЗАПИСКИ ЦАГИ 20 0 7

№ 1 — 2

УДК 533.6.013.2.011.3/5:629.7.025.73

НЕСТАЦИОНАРНЫЕ АЭРОДИНАМИЧЕСКИЕ НАГРУЗКИ НА ПРОФИЛЕ В ПРОИЗВОЛЬНОМ ВЕРТИКАЛЬНОМ ПОРЫВЕ И ИХ МОДЕЛИРОВАНИЕ С ПОМОЩЬЮ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

А. Н. ХРАБРОВ

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

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

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

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

[5],

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

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

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

1. Исследуется движение двумерного тонкого симметричного профиля (пластинки) с постоянной скоростью в идеальной несжимаемой жидкости. Задача решается в безразмерном виде. При этом предполагается, что скорость движения профиля в жидкости и = 1 и длина его хорды Ь = 1. Решение рассматривается в связанной с профилем системе координат, начало которой расположено в точке, совпадающей с его условным центром тяжести. Относительно этой же точки рассчитывается аэродинамический момент тангажа. Профиль в этой системе координат расположен вдоль оси х от точки х0 до точки х0 - 1 (х0 — центровка профиля).

Нормальная скорость на профиле в случае его проникновения в ступенчатый порыв имеет вид:

Г-1, X е [х0 -1,х0],

V« (X, 0 = ^ ’ 10 ’ 0Ь (1)

[0, х е [х0 -1, х0 - (].

При этом предполагается, что в начальный момент времени носик профиля начинает проникать в вертикальный порыв.

Аналогично работе [6] на основании решения соответствующей краевой задачи для комплексно сопряженной возмущенной скорости получим:

ём> 1

“г п ^

г - хп +1

г - хп

(2)

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

1;2(х)&гл=- т-х-ч ^ (3)

которое может рассматриваться как интегральное уравнение для функции и2, так как в правой части стоит известное выражение, зависящее от распределения вертикальной скорости в порыве. Это интегральное уравнение Вольтерра первого рода с ядром, содержащим интегрируемую сингулярность. В работе [5] на основании алгоритма [7] была разработана и проверена эффективная методика решения таких интегральных уравнений. После нахождения функции и2 нестационарные нагрузки, действующие на профиль, могут быть рассчитаны с помощью формул Блазиуса — Чаплыгина. В данном случае получим для коэффициента подъемной силы:

^-'4“ х0 ^(х0 -0)(1 -х0 + 0) “0 + 21 1( и,2(0) • (4)

Мх0-1 0>/(х0 -1 -0)(х0 -0)

Аналогично для коэффициента аэродинамического момента тангажа

*0 I _ Е , -0 __________________

mz = _2 f Vn (1 + 2 Ш -Г-0--7 d Е_2 “Г f Vn W( x0 _Е)(1 _ x0 + Е) d E_

*0 _i V1 _-0 +E dt-0 _i (5)

_2 ( -0 _ Щ ^( -, _E)(1 _ -0 +Е) dE + 2 ( -0 _ i) jv( -0 _“i2^0 _е) .

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

S = Су1 + СУ2 + Су3, (6)

mz = mZ1 + mz 2 + mz 3,

где члены Cy1 и mZ1 представляют собой квазистационарные составляющие. Через Су2 и mZ2 обозначены члены, содержащие явные производные от времени от интегралов по профилю, которые по аналогии с движением профиля как твердого тела можно назвать присоединенными массами. И наконец, члены cy3 и mz3 выражаются через интегралы по вихревому следу. Они зависят от предыстории движения и описывают нестационарные эффекты запаздывания.

2. В случае ступенчатого вертикального порыва подстановка выражения для нормальной скорости на профиле (1) в правую часть интегрального уравнения (3) после вычисления определенных интегралов приводит к следующему выражению:

f (t ) =

— + —arcsin(2t _ 1) _1J t (1 _ t ), t e[0,1],

4 2п n (7)

П ri -1

2, t e [1, œ]

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

су1 = п + 2агс8т(2^ -1) - 4^Jt(T—tj.

Производная по времени от второго интеграла для случая проникновения в единичный ступенчатый порыв (присоединенная масса) имеет следующий вид:

Су 2 = 4л/?0'-0.

Таким образом, полный нестационарный коэффициент подъемной силы определяется формулой:

су = п + 2агс8т(2^-1) + 2Г .—и2(^)“^ =. (8)

у 0

Аналогичные вычисления для коэффициента момента тангажа дают в результате выражение:

*'=*(х0 - 4)+2 (х0 - 4 ) агс5‘"(2'-1)+2(х0 - 4) ^ -,)

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

Рис. 1. Результаты решения задачи о постепенном входе профиля в вертикальный порыв ступенчатого вида:

1 — Су1 + Су2, 2 — Суз; 3 — полный су, 4 — точное аналитическое решение

Кюсснера

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

Результаты вычисления составляющих коэффициента подъемной силы по результатам численного решения интегрального уравнения Вольтерра представлены на рис. 1. Составляющие Су1 + Су2, зависящая от развития вихревого следа часть су3 и полное значение коэффициента подъемной силы Су, обусловленное проникновением профиля в ступенчатый порыв, показаны различными типами линий. Там же показана функция Кюсснера, представляющая собой точное аналитическое решение данной задачи. Эта функция выражается в виде интеграла от Бесселевых функций и построена по ее табличным значениям, приведенным в [3]. Видно, что полученное численное решение прекрасно совпадает с классическим точным решением. Исключение составляет окрестность точки ^ =1, где в численном решении возникает некоторая неустойчивость, выражающаяся в появлении высокочастотных колебаний. Это обусловлено особенностью в точке ^ = 1, которая существует у функции (7), определяющей правую часть интегрального уравнения Вольтерра. Возникновение такой колебательной неустойчивости является характерным свойством некорректных задач, которыми являются интегральные уравнения первого рода.

При вычислении интеграла по вихревому следу в выражении для коэффициента Су использовались квадратурные формулы разложения в окрестности сингулярности, полученные в работе [5].

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

Г ) = ^~Гd Х = \~Г -Т^ Т, (9)

а т J d т

о о

где 0 — ступенчатая функция

0(ґ) =

[0, t є (-да, 0),

І1, t є [0, да).

В линейной задаче проникновения профиля в произвольный вертикальный порыв имеем:

Здесь ^) — профиль вертикальной скорости в порыве; переходная функция I ^) —

реакция на элементарный ступенчатый порыв (функция Кюсснера). Если аппроксимировать функцию Кюсснера рядами экспонент и оставить только два члена, то

1с (0 = 2пГ 1 -Ье~в1 -(1 -Ь)е~в*

(11)

Здесь Ь, Р1 и Р2 — некоторые постоянные, которые в дальнейшем будут идентифицированы. При этом было учтено, что 1с (0) = 0. Подстановка (11) в (10) дает в результате выражение:

с ^) = 2пуг ^) - 2пЬ \——е в1<^ т) а т- 2п(1 - Ь)\——е в2^ т) а т. у — J а т і а т

0 0 С использованием формулы интегрирования по частям можно записать также:

г г

су (г) = 2 лЬр11 у^е_в1(г-х) й т - 2 п(1 -Ь)Р21 Уёе~^2('1~т') й т.

(12)

(13)

1

После преобразования Лапласа получим:

Су (р) = 2пЬР1— (р) Р + 2п(1 - Ь)$2— vg (р)

У Р g Р + Р1 Р g Р + Р2

Это выражение во временной области эквивалентно динамической системе:

X +Р1-Х1 = 2пЬвlVg (г),

X2 +Р2X2 = 2п(1 -Ь)P2Vg (г),

Су = Х1 + X2.

(14)

Здесь Х1 и Х2 — некоторые внутренние динамические переменные. Аналогичная динамическая система с другими правыми частями использовалась ранее для моделирования нестационарных нагрузок на профиле при изменении угла атаки и угловой скорости тангажа [5], а также угла отклонения элерона [6]. В связи с этим константы р1 и р2 резонно выбрать, соответствующими полученным ранее значениям, основанным на аппроксимации функции Теодорсена. При этом ввиду линейности задачи возможно сохранить одну и ту же динамическую систему для моделирования всех нестационарных эффектов, объединив слагаемые в правой части. Таким образом, для аппроксимации функции Кюсснера необходимо идентифицировать только параметр Ь. Это можно сделать с помощью минимизации функционала:

ф=!

(15)

с,

4

б

2

О

5

10

15

Рис. 2. Аппроксимация решения Кюсснера динамической системой второго порядка:

1 — Х\, 2 — Х2; 3 — полный Су = Х\ + Х2; 4 — решение Кюсснера

В этом выражении для функции Кюсснера введено обозначение K(t). Были проведены вычисления, в которых в качестве tj брались 200 точек равномерно распределенных в интервале t е [0, 20]. Для значений Pi = 0.1673 и р2 = 0.8598 [5] при минимизации с помощью симплексного метода Нелдера — Мида получено значение b = 0.2710 с конечным значением функционала Фтт = 0.1219. На рис. 2 показаны полученные результаты для аппроксимации функции Кюсснера. Разными типами линий показаны переходные процессы для Xi(t) и X2(t), их суммы Cy(t), а также точное значение функции Кюсснера. Видно, что результаты аппроксимации вполне удовлетворительны.

Если данная точность аппроксимации в какой-то задаче недостаточна, следует повысить порядок динамической системы

Значениям Р1 = 0.07237, р2 = 0.37980 и Рз = 1.36407 [5], полученным при аппроксимации функции Теодорсена для величин В1 и В2 в результате минимизации (15), соответствуют значения В1 = 0.1394 и В2 = 0.2687. В этом случае конечное значение функционала равно Фтт = 0.0305, что существенно меньше полученного минимального рассогласования для второго порядка.

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

ICy (t) = 2п[і - B1e~eit - B^2 - (1 - B - B2)e~e3t ,

(16)

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

Xi +Pi Xi = 2nBiPiVg (t),

X2 +P2X2 = 2nB2?2vg (t)

X3 + РэX3 = 2n(i - Bi - B2 )Рзvg (t), Cy = Xi + X2 + X3.

(i7)

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

(х) = <

0, х е (-да, -Ь/2], 0.5 + х/Ь, х е[-Ь/ 2, Ь/ 2],

1, х е [Ь/2, да).

(18)

Для второго вида порыва использовалось следующее выражение:

0, х е (-да, -Ь/2],

0.5 + 0.5оо8

2пх Ь ’

х е[-Ь/ 2, Ь/ 2],

1, х е [Ь/2, да).

Рис. 5. Различные составляющие точного численного решения

Рдафгршштнйнтсрявяоео урав&еябо ШлееягёгдакефуйВДВДрзрйЬа

Рис. 4. Различные составляющие точного численного решения для градиентной ступеньки Ь = 1:

1 — Су 1; 2 — Су2 3 — су3; 4 — полный су, 5 — аппроксимация с помощью динамической системы второго порядка

В обоих случаях Ь — характерная длина порывов. Вычисления были проведены для значений Ь = 1 и Ь = 3. На рис. 3 показаны результаты решения интегрального уравнения (3) для

М2

в случае профиля, проходящего через порывы (18) и (19) различной длины Ь. В момент времени ^ = 0 носик профиля во всех случаях проходит точку х = -Ь/2. При решении интегрального уравнения для порывов заданной произвольной формы правые части (3) находились численно. Также численно проводилось и интегрирование нестационарных аэродинамических нагрузок в соответствии с выражениями (4) и (5).

-2

О 2 4 /6

Рис. 7. Различные составляющие точного численного решения для одного периода косинуса L = 3 (обозначения те же, что на рис. 4)

На рис. 4 показаны полученные составляющие для коэффициента подъемной силы при прохождении порыва в виде градиентной ступеньки длиной L = 1. Различные типы линий соответствуют составляющим точного численного решения (Cyi, Cy2, Суз) и суммарному значению

Cy(t).

Там же представлены результаты приближенного математического моделирования Cy(t) с использованием динамической системы второго порядка (14). Видно, что простая математическая модель удовлетворительно описывает результаты, полученные на основе численного решения интегрального уравнения Вольтерра. Данные на рис. 5 соответствуют аналогичным результатам сравнения для порыва ветра той же формы длиной L = 3. Видно, что по мере увеличения длины порыва точность аппроксимации увеличивается.

Рис. 6 и 7 представляют аналогичные результаты сравнения для ветрового порыва (19) длиной L = 1 и L = 3. Видно, что и в этом случае математическое моделирование с помощью динамической системы второго порядка дает результаты для развития нестационарного коэффициента подъемной силы в качественном и количественном отношении хорошо согласующиеся с суммарным Cy(t), полученным на основе точного численного решения.

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

Работа выполнена при финансовом участии РФФИ (проект 03-01-00918), а также Президентского гранта поддержки ведущих научных школ НШ-2001.2003.1 (руководитель академик Г. С. Бюшгенс).

ЛИТЕРАТУРА

1. Küssner H. G. Zusammenfassender Bericht über den unstationären Auftrieb von Flügeln // Luftfahrt-Forshung. — 1936. Vol. XIII, N 12.

2. Sears W. R. Some aspects of non-stationary airfoil theory and its practical application //

J. of the Aeronautical Sciences. — 1941. N 8.

3. Некрасов А. И. Теория крыла в нестационарном потоке. — М. — Л.: Изд.

АН СССР. — 1947.

4. Седов Л. И. Плоские задачи гидродинамики и аэродинамики / Изд. 3-е, переработанное. — М.: Наука. — 1980.

5. Храбров А. Н. Математическое моделирование влияния схода вихрей на нестационарные аэродинамические характеристики профиля при его произвольном движении // Ученые записки ЦАГИ. — 2002. Т. XXXIII, № 3 — 4.

6. Храбров А. Н. Математическое моделирование нестационарных аэродинамических нагрузок на профиле при динамическом отклонении органа управления на задней кромке // Изв. РАН. МЖГ. — 2005, № 3.

7. Hairer E., Lubich Ch., Schlichte M. Fast numerical solution of weakly singular Volterra integral equations // J. of Computational and Applied Mathematics. — 1988. Vol. 23.

Рукопись поступила 31/12006 г.