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

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

Том XXXIV

20 0 3

№ 1—2

УДК 629.735.33.015.073

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

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

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

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

Исходными данными для задания случайных скоростей порывов являются форма спектральной плотности и два параметра — масштаб и дисперсия. Наиболее распространенными являются модели спектральной плотности Драйдена и Кармана.

В случае изотропной турбулентности энергетические функции спектральной плотности для этих моделей имеют вид [1]: модель Драйдена

В. П. КУЗЬМИН

\¥ = (и, V, ч>)Т = W(х, у, z, X).

1+(ьо)2 ]3

модель Кармана

где Ь — масштаб; с — среднее квадратическое отклонение; 0 = (0х, Оу, О2) — вектор пространственной частоты, О2 = О^ + О2 + О2, ак = 1,339 — постоянная Кармана.

т

Функция Е(О) позволяет определять спектральные плотности компонентов случайного вектора скорости порыва. Трехмерная спектральная плотность для /-го компонента вектора скорости ветра определяется по формуле [1]:

5 (о о а)-£<П)-о2)

Sj X , аy , аz) -

4пО4

Спектральные плотности от меньшего числа переменных определяются путем интегрирования:

Sj (о X, о z )- j Sj (о X, о y, о z ) d о y ,

(О х )= | (О х, О г ) й О г .

—да

Конкретные формулы для различных моделей турбулентности приведены в [12]. В приземном слое атмосферы (к < 300 м) параметры спектральных плотностей различных компонентов скорости ветра различны и изменяются с изменением высоты [2]. Использованная в данной работе модель взята из [3]:

СТ = Гл т л ою/1 /-,лл\~|0,4 ' ^ = к , ^и = ^ = ^ (аы/V) .

[0,177 + 0,823 (h¡300)]

Будем считать ортогональные компоненты порывов некоррелированными. Тогда известные спектральные плотности для различных компонентов позволяют задавать случайные реализации скоростей порывов. В данной работе учитывается распределение скоростей порывов по длине (х) самолета и размаху (z). Распределение скоростей порывов по высоте самолета не учитывается, поскольку считается, что высота самолета значительно меньше его длины и размаха крыла.

Наиболее общим способом задания реализаций является представление скоростей ветра в виде рядов Фурье. В случае двумерной турбулентности компоненты скорости порыва могут быть заданы в виде

M N

W<X, z) - ^^ Alk cos (о*х + Oz/z) +Blk sin (о„-х + о^), (1)

j=1 i-1

где A/j и B/j — независимые случайные гауссовские величины с дисперсиями

D[A/j ] - D[By ] - Sj (ох/, оzj) ДохДоz , к = 1, 2, 3.

Для модели Драйдена одномерные реализации скоростей порывов могут задаваться с помощью стохастических дифференциальных уравнений (формирующих фильтров) [11].

2. Определение нестационарных аэродинамических характеристик самолета. Для

анализа воздействия турбулентных порывов на самолет существенным является распределение скоростей порывов по самолету. Для определения нестационарных аэродинамических характеристик используется панельный метод [5] — [9]. Летательный аппарат представляется системой плоских несущих поверхностей, параллельных вектору скорости набегающего потока. Рассматриваются малые гармонические колебания несущих поверхностей с частотой ю. Нормальная к несущей поверхности безразмерная скорость потока (скос) wn связана с амплитудой колебаний несущих поверхностей

W д f

wn ---f + ikfn, (2)

V0 д x

со

где /п = (/ ■ п) скалярное произведение вектора перемещений / = (/у, /2) и вектора нормали п = (пу, п2) ; ¥0 — скорость набегающего потока; к = ю/2/¥0 — приведенная частота;

¡2 = Ц 2 — полуразмах крыла, к которому отнесены также все линейные размеры самолета. Основное интегральное уравнение для определения давлений имеет вид [5] — [7]:

(х0. Уo, 20 ) = ПЦ К (х0. Уo, 2о. П. С. к ,M )Ар. п. сМ . 5 (3)

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

N

^ = Е Кч Аср . (4)

у = 1

где

Ку =4 ДО К (^ .Пу .С/у,к ,м) *у

5

= х0/ - Ху. Цу = у0/ - у у. ^у = 20/ - 2 у — переменные координаты точек на у-ой панели,

отсчитываемые от -ой контрольной точки и взятые с обратным знаком.

Вектор скосов в контрольных точках всех панелей определяется по формуле (2) для заданной формы колебаний, после чего вызванные этими колебаниями безразмерные разности давлений на панелях определяются из решения системы линейных алгебраических уравнений (4):

Аср = К "Ч. (5)

где К = |К ||.

Безразмерные коэффициенты обобщенных сил определяются путем интегрирования по всем несущим поверхностям

е =—1Г /

^пт „ II У к

„ т Асрп^-5о

где /т — перемещение несущих поверхностей по т-ой форме; Асрп — разность давлений,

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

Для определения стандартных продольных и боковых аэродинамических характеристик расчеты проводились при числе Маха, равном 0,2, и малом значении приведенной частоты к □ 1 путем задания соответствующих форм колебаний.

Расчеты аэродинамических коэффициентов и приращений сил и моментов от порывов ветра проводились для тяжелого самолета типа Ил-96 с размахом крыла I = 58 м и площадью крыла

5о = 350м 2.

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

Л={0}'Л={0} • Безразмерные аэродинамические коэффициенты определяются соотношениями:

сау = Яе (вп), та = Яе (022),

с« = Яе|в1! 1, та = Яе Г вх

а , т

су + су7 = 1т

к в12

2 ' 2

а , т

. ,, т7 + т77 = 1т|

у у 1 к У 2 2 I к

в22

В боковом канале задаются три антисимметричные формы:

={0{' ={—,{• * = {—

Безразмерные аэродинамические коэффициенты определяются соотношениями:

св =

= Яе (в12), т\ = Яе (в22), тв = Яе (вз2),

в =

св=— Яе | % 1, тв =— Яе Iв^ 1, тв = — Яе ( % |,

св+ с» = 1т Iе12 1, т в+ т Г = 1т

к Г

у "У

в221, = 1т Г вз2

с»* = 1т в |, тг* = 1т | ^ |, т»* = 1т Г вз3

к у у у к у х У к

В табл. 1, 2 приведены продольные и боковые аэродинамические характеристики.

Таблица 1

Продольные аэродинамические характеристики

са Г7 . а с 7 + с а с г7 . а — с 7 + с

сУ 5,90 1,68 -0,07 -1,85

тг -0,46 -1,28 -0,30 0,68

Таблица 2

Боковые аэродинамические характеристики

св

Гу , в с у + с4 сг* св с»У — св

с7 -0,75 -1,04 -0,14 -0,37 -0,30

тх -0,25 -0,12 -1,03 0,005 -0,62

ту -0,25 -0,58 -0,02 0,02 -0,13

3. Определение приращений сил и моментов от порывов ветра. Для определения

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

атмосферы, достаточно определить передаточную функцию от порывов ветра к безразмерным аэродинамическим коэффициентам.

Рассмотрим эту процедуру для возмущенных скоростей, нормальных к несущим поверхностям, т. е. для боковых и вертикальных порывов. В данном случае передаточная функция зависит от двух безразмерных частот, kx = Qи kz = Qz/2, и гармонический скос единичной амплитуды имеет вид

wn (x, z, t) = wn (x, z)e'kxt = -e'kxt ex pikxx + ikzz) .

Откуда получим, что для поперечных составляющих порывов скос задается в виде

wn (x, z) = - exp (ikxx + ikzz) . (6)

Коэффициенты аэродинамических сил и моментов будут иметь вид c = H (kx, kz ) e'kxt, где H (kx, kz ) = H1 (kx, kz ) + iH2 (kx, kz ) — комплексная передаточная функция.

Для учета влияния продольных порывов считается, что выполняется интегральное соотношение (3) для размерных переменных. Тогда при изменении продольной скорости получим

N

aw, (j - U ) = wn, = ^ KJACPJ (J + UJ ), J=1

и изменение коэффициента давления, обусловленное продольными порывами, в линейном приближении определяется соотношением

N

Acp, = (j - u, ) 2 K-1 (kx ) awj (j -UJ ) - Acp0, , j=1

* dfn где Acp0 = K aw — стационарное распределение давления, а aw =--- — углы наклона

д x

несущих поверхностей к набегающему потоку.

Очевидно, что если порыв стационарный (u = const), то изменение любого аэродинамического коэффициента определяется формулой Acpu = -2Acp0u , соответствующей

линейному приращению аэродинамических сил при изменении скоростного напора.

Отметим, что в данной работе размерные обобщенные аэродинамические силы

oV 2

представляются в виде Q = c ^ So, поэтому всякое изменение обобщенной силы трактуется как

изменение аэродинамического коэффициента Ac = 2 AQ .

pVo2So

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

u = - exp (ikxx + ikzz) .

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

Известные передаточные функции позволяют определять спектральные плотности сил и моментов:

Sc (kx, kz ) = Re \h (kx, kz) HT (kx, kz)] Sw (kx, kz),

Sc (к* ) = { Зс (к*, К) йк2

Рис. 1

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

М N

Щх) = ^^ |Н/1 [ А сое (кХ1х + фу.) + Вг] эш (кХ1х + фу.)] :

/ = 1 г=1

где

ф/ = ЗГ^ [Н2 (кхг, к2/ )/Н1 (кхг, К/ )] •

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

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

4. Учет конечных размеров самолета в первом приближении. Для удобства практического использования метод приближенного учета конечных размеров самолета должен удовлетворять двум основным условиям:

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

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

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

= 0 + ™пХ + wnz ,

а аэродинамические коэффициенты определяются формулой

(7)

С = С п^п0 +

д с

д7х

дс д ^

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

д с д с

порывов и определению соответствующих аэродинамических коэффициентов с^п,-,- для

дwx дwz

п п

трех компонентов скорости порыва и пяти аэродинамических коэффициентов.

4.1. Определение коэффициентов влияния. Рассмотрим вначале задачу определения коэффициентов влияния порывов ветра и их градиентов на изменение аэродинамических сил и моментов, т. е. вторую из перечисленных выше задач. В табл. 3 приведены безразмерные коэффициенты влияния различных компонентов порывов ветра и их производных на приращения аэродинамических коэффициентов, определенные путем расчета передаточных функций для малых значений частоты кХ = 0,05.

Таблица 3

Коэффициенты влияния порывов ветра

си дс дих дс д^ с" дс дVх дс дV дс дwx дс дwz

СУ -3,02 0,41 0 0 0 -0,21 5,90 -1,86 0

Сz 0 0 0,08 0,75 -0,30 0 0 0 -0,04

тх 0 0 0,64 0,25 -0,13 0 0 0 -1,01

Шу 0 0 0,12 0,25 -0,62 0 0 0 0,06

mz 0 0,04 0 0 0 0,03 -0,46 0,69 0

Как следует из анализа результатов расчетов передаточных функций, достаточно воспроизвести 8 случайных величин и0, и1, "0, Vх, V1, w0, wx, wz и определить зависимости

Су (и0, Wо, wx), Сz (V, Vх), Шх (V), Vх, иг, wz), ту (V, Vх), тг (Wо, ).

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

Безразмерные коэффициенты давления определяются соотношением (5). Если частота мала,

то коэффициент давления в первом приближении определяется по формуле:

Лср «М + 7к^2 ) . (8)

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

При малых частотах скос от вертикальных порывов представляется в виде

wn = w «-(1 + ¡кхХ + 7кг2) .

Местные давления в первом приближении определяются соотношением

Лср = - (М1 + 7кхМ2 )(1 + 7кхх + 7кг2) = -М1 - 7 (кхМ2 + к^х + кгМ^) . (9)

При повороте самолета по тангажу вокруг оси 02 скос и коэффициенты давлений определяются соотношениями:

w « -1 + 7кхх ,

Лср =-М1 - 7кх (М2 - М1 х) . (10)

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

Сравнивая соотношения (9) и (10), легко определить, что вертикальный ветер эквивалентен углу атаки, а градиент вертикального порыва по продольной оси эквивалентен угловой скорости -ю2 + а и, следовательно, коэффициент пропорциональности между вертикальным ветром и продольным градиентом вертикального ветра и каким-либо аэродинамическим коэффициентом определяется через соответствующие аэродинамические коэффициенты:

— =са, (11)

дw

— = -с5г + са. (12)

дwx

Вращение самолета вокруг оси ОХ создает скосы: на горизонтальных поверхностях — w = wn = 7гюх, на вертикальных поверхностях — /= wn = 7уах , тогда как вертикальные порывы создают скосы только на горизонтальных поверхностях. Если влияние скосов на вертикальных поверхностях на изменение аэродинамического коэффициента мало, то можно принять

дС ^ . (13)

д^

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

Наиболее сложным является установление соответствия между угловой скоростью рысканья и градиентами порывов [4]. Вращение самолета вокруг оси ОY создает скосы,

аналогичные боковому ветру и двум градиентам V и и2 .

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

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

дс dv

= Л

дс

ä7

=

(с, = 0)- сг

(14)

(15)

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

дс = сЕу (су )-с5у (су = 0). Влияние постоянного продольного ветра определяется очевидным соотношением

(16)

дс

— = -2с . ди

(17)

Таким образом, соотношения (11) — (17) приближенно определяют расчетные коэффициенты влияния, значения которых приведены в табл. 4. В ней выделены коэффициенты влияния, значения которых отличаются от точных значений, приведенных в табл. 3.

Таблица 4

Расчетные коэффициенты влияния порывов ветра

си дс дс duz с" дс дс д" с" дс дс д^

сУ -3,02 0 0 0 0 0 5,90 -1,86 0

сг 0 0 0,08 0,75 -0,30 0 0 0 - 0,14

Шх 0 0 0,64 0,25 -0,13 0 0 0 -1,03

Шу 0 0 0,12 0,25 -0,62 0 0 0 - 0,02

mz 0 0 0 0 0 0 -0,46 0,69 0

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

4.2. Схема вычисления порывов ветра и их производных. Если представить порывы ветра точно в виде (7), то дисперсии производных порывов по оси ОХ и зависящих от этих производных аэродинамических коэффициентов, при использовании моделей Драйдена или Кармана, будут бесконечными. Для того чтобы избежать недифференцируемости случайных процессов, определяющих производные порывов, в работе [4] предложено вычислять градиенты скоростей порывов как конечные разности для некоторых характерных точек на самолете. Точка 0 расположена в центре масс самолета, точки 1 и 2 — симметрично на обеих половинах крыла с расстоянием между ними, равном 0,71, точка 3 — в центре давлений на горизонтальном

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

Производные скоростей порывов по переменным x и z определяются соотношениями:

Wz = (W2 - W1 )/2Az = AW/2Az , Wx = (W0 - W3 )/Ь ,

где W3 = H3 (kx )Wo, Ь = Ьг для вертикальных порывов и Ь = Ьв для боковых порывов.

Точное значение передаточной функции H3 для замороженной турбулентности имеет вид

Hз (kx ) = e~kxb , (18)

а спектральная плотность для разности AW определяется интегралом

то

Saw (kx) = 2 j(l - cos (kz2Az))Sw (kx, kz)dkz . (19)

0

Величины W) и AW являются некоррелированными.

В работе [12] приведены выражения для спектральных плотностей разностей порывов в двух точках для различных компонентов порывов и различных моделей турбулентности. Эти выражения, как и передаточная функция (18), не являются дробно-рациональными и, следовательно, не представляются с помощью дифференциальных уравнений.

Использование конечных разностей для вычисления производных порывов оказывается недостаточным для удовлетворительной аппроксимации спектральных плотностей сил и моментов в широком диапазоне частот для приземного слоя атмосферы, где масштабы турбулентных порывов становятся соизмеримы либо меньше характерных размеров самолета. Это означает, что в данном случае нельзя ограничиться только аппроксимацией выражений (18) и (19) с помощью дробно-рациональных выражений, а необходимо (для малых масштабов порывов) использовать дополнительно различные эмпирические зависимости. Эти зависимости подбираются непосредственно из условия обеспечения приемлемой точности аппроксимации соответствующих спектральных плотностей для безразмерной частоты k < 3, что достаточно для моделирования движения неманевренного самолета как твердого тела на режиме посадки.

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

Формирующие фильтры или спектральные плотности для величин w и wx подбираются из условия совпадения спектральных плотностей подъемной силы и продольного момента, а спектральная плотность для производной wz — путем аппроксимации спектральной плотности момента крена.

В результате подбора были найдены соотношения:

H з (kx ) =

1 + kx>Ь +1/2 (kx>Ь )2

Sw0 (kx ) = SW (kx ) :

о (L )_ 2 aAwLAw SAw \kx)~'

п 1 + (kxLAw )2

(20) (21) (22)

1

Ь. 2 2с.

ГДе ^ = [1 + 0,32^/Аг]' СД. = 1 + 0,32Ь./Аг '

Приращение угла атаки и производные по переменным х и г определяются соотношениями:

Аа = .оР (Аг/Ь.) , = .х = (Аг/Ь. ), ^ = = (Аг/Ь.) , (23)

эмпирические множители в соотношениях (23) учитывают эффект осреднения распределения порывов по размаху крыла при малых масштабах и имеют вид

Ра (ЧЬм>) = ^1 - ехР (-Ь.12Аг), Р ( ) = 11 - ехр (-Ь.1 Аг),

Г1 при Ь. >, Аг (24)

Р (Аг/Ь. ) = ] ^^

и Ь.1 Аг при Ь. < Аг.

Из приведенных выше соотношений легко заметить, что величина Аг может быть выбрана произвольной, поскольку ее изменение приведет лишь к изменению эмпирических множителей. Принятое в данной работе значение Аг = 0,35/ отличается от других возможных значений тем, что в этом случае при больших значениях масштаба эмпирические множители (24) становятся близкими к единице, а спектральная плотность для производной близка к спектральной плотности для разности порывов (19).

Продольные поршЬ»

Рис. 2

10S

Рис. 3

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

5и ( кх ) =

я2 ,-т2;

р и = ри <

2 4 р2 / ~ \2 ' П 1 + (^иК )

0,44 + 0,56ехр

^и = ии

2 - ехр

( А^ >

V ии ;

-1,1

- 2 — и

(25)

5

Аи (кх) =

г _2

^Аи° Аи

П 1 + (ЬАикх )

ЬАи = Ьи

1 - ехр

о лАг

-3'4 т

и у

Аи

= 2°2 I1 - ехР

"2,8 Г

. (26)

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

Вертикальные порыЬь

БокоЬне порыЬы

Рис. 4

Сравнение точных и приближенных спектральных плотностей для основных коэффициентов показано на рис. 2—5. Для удобства сравнения при различных частотах на

Боковые порыЬы

Рис. 5

рисунках приведены амплитудные характеристики (V 5 ), причем точные значения — сплошными

линиями, а приближенные — штриховыми.

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

и0 = -и0 (У/ии ) + ри\¥1ии 1и ,

эДи :

М = -Дыг (V/LM ) + Сд„ {ЩГДи^д w о = -( w0 +zw )(V/LW ) + aw¡3^L~^w, = -zw (V/Lw ) + aw (l -1//3 У 3V/Lw ^

^z = -Дwz (Vl^w ) + ^/^¿^Ди , w2 = -(w2 - xw)(V/b ), xw = -(w2 + Xw - 2w0)(Vlb ),

v = - (v + z v) (V/Lv) + a v 3V/L¿v,

¿v = -Zv(V/LvV + a v(l - 1/^3 3V¡L£v,

V2 = -(V2 - xv)(V/b ), Xv = -(v2 + Xv-2vo)(V/b ),

ДVz = -Д v¿ (V/Lv) + a WL^ v.

В приведенных уравнениях переменной Е, с различными индексами обозначены независимые гауссовские белые шумы единичной интенсивности, а все параметры в этих уравнениях определяются формулами (22) - (26).

Из данных переменных формируются 7 случайных величин:

w = w0Fa (Lw ) C0s Y + v0Fa (L v) sin Y, v= vcos Y + w0sin y,

и = u00,

= ДwFx (Lw ) cos y + Д vFx (L v) sin y ^ 2Дг '

f = ( w2 - w0 ) F¿ (Lw ) cos Y + ( v - v0 ) F¿ (Lv) sin Y ra¿g h '

<» ^gi =

_ У - v ) cos Y + ( w2 -wp ) sin Y

Ьв

yg 2 = -

Ди

где Fa(L) = J 1 - exp (-L/2Дг) , F¿ (L) = J 1 - exp (-L/Дг) , Fx (L) = ^

2Аг

Г 1 при X > Аг при X < Аг.

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

формулам:

Сх = Сх (их -u, иу - V и1 -w, а + Юzg, Р - Юyg1, Юх + Юxg, Юу, -Юzg ), Су = Су (их - u, иу иг -w, а + Ю^, Р - Юу^ Шх + Юxg, Юу, Ю z -Юzg ),

Сг = Сг (их - u, иу - V иг -w, а + у, Р - Юх + Юxg, Юу, Ю z - Юzg ) +

+ ^у (Су = 0)ЮУ51 + у - ^у (Су = 0)]Юyg2 ,

тх = тх (их-и, иу - V, и2-w, а + Юzg, Р-Юу81, Юх ,

+ т1У (Су = 0)юygl + туу - туу (Су = 0)^

Юу, Юz -Юzg) +

|Ю Уg 2 ;

ту = ту (их-и, иу - V, Uz-w, а + ю^, р-ю„„,, Ш+Ю, ю„, Ш-Ю ) +

JygЪшx^шxg'шy'шz zg ,

■ту'

(Су = 0)Юу51 + туу - туу (Су = 0)]Юyg2 :

т2 = mz (их - и, и у - V, Uz - w, а + у zg, Р - ю ygl, ю х + ю xg, ю у, ю z - ю zg ) •

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

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

6. Оценки точности движения самолета при автоматическом управлении. Для оценки влияния точности аппроксимации спектральных характеристик сил и моментов от порывов ветра на статистические характеристики движения самолета рассматривается линейная модель автоматической посадки самолета типа Ил-96 с постоянной скоростью.

В боковом канале алгоритм управления самолетом сводится к стабилизации заданной линии пути, а в продольном канале — к стабилизации номинальной траектории в вертикальной плоскости, которая соответствует движению по прямолинейной глиссаде с последующим экспоненциальным выравниванием [10].

Расчеты проводились для скорости V = 70 м/с и интенсивности турбулентности, соответствующей aw = 1 м/с.

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

Расчет одной реализации производится для одного компонента порывов и для заданных значений частот Ох1 и О^ . Отклонения аэродинамических коэффициентов задаются в виде

гармоник

ас = в'Ох'хИк (О х, О д ,

где к =1, 2, 3, / =1, 2, ..., N у = 1, 2,...,М.

Расчет передаточных функций проводится на сетке N *М = 40 * 50 в диапазонах Ох/2 е [0,4] и О¿2 е [0,8].

о

■8

а

1*Г

гЬа

"Т|о гЬо

гЬо

N

3

о

Яэ'

_ __

1Г

"ТЙо 1^0"

"гЬо

Ьысота ( м )

Рис. 6

-Зг

"гЬо

Ьысота ( м )

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

3 N М

куу=ХХХ^/З-

к = 1 I = 1 у = 1

Для приближенной оценки средних квадратических отклонений расчет проводится по методу Монте-Карло для 1000 реализаций. При этом порывы ветра и их производные задаются с помощью формирующих фильтров, а линейные приращения аэродинамических коэффициентов определяются с помощью коэффициентов влияния, приведенных в табл. 4.

Сравнение точных и приближенных значений средних квадратических отклонений различных фазовых координат приведено на рис. 6 и в табл. 5.

Таблица 5

Сравнение точных и приближенных средних квадратических отклонений

Точные значения Метод Монте-Карло

И, м 200 90 30 5 200 90 30 5

Пу 0,056 0,061 0,051 0,039 0,060 0,058 0,052 0,039

ю z, град./с 0,67 0,66 0,48 0,31 0,68 0,67 0,50 0,30

град. 1,14 1,04 0,78 0,51 1,17 1,04 0,78 0,49

Уу, м/с 0,73 0,74 0,61 0,46 0,72 0,76 0,61 0,45

АН, м 1,65 1,75 1,69 1,32 1,68 1,67 1,61 1,22

ю х, град./с 0,63 0,81 1,04 1,15 0,66 0,80 1,05 1,09

у, град. 0,83 1,03 1,25 1,31 0,84 0,99 1,22 1,26

У, м/с 0,38 0,46 0,54 0,50 0,36 0,45 0,51 0,50

АХ, м 1,84 2,15 2,41 2,38 1,65 1,94 2,12 2,16

ю y, град./с 0,29 0,36 0,44 0,41 I 0,28 0,34 0,43

Y, град. 0,83 0,99 1,08 0,85 | 0,78 0,92 1,03

0,42 0,83

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

Отметим, что эти ошибки не связаны полностью с использованием метода Монте-Карло, поскольку для объема выборки в 1000 реализаций относительная величина среднего квадратического отклонения оценок составляет « 2,2% и, следовательно, относительные ошибки более

5—6% являются значимыми. Однако ошибки в 10—15% при моделировании движения самолета вблизи земли следует считать небольшими, поскольку любая из моделей турбулентности в приземном слое является приближенной.

Работа выполнена при поддержке Российского фонда фундаментальных исследований, проект(01-01-00431).

ЛИТЕРАТУРА

1. E t k i n B. Dynamics of atmospheric flight//John Wiley & Sons.— Inc. 1972.

2. Характеристики ветровых возмущений в нижних слоях атмосферы//Обзор ОНТИ ЦАГИ.— 1979, № 545.

3. Shakarian A. Application of Monte-Carlo techniques to the 757/767 autoland dispersion analysis by simulation//AIAA 83-2193.— 1983.

4. E t k i n B. Turbulent wind and its effect on flight//J. of Aircraft.— 1981. Vol. 18, N 5.

5. Rodden W. P., Giesing J. P. and K a l m a n T. P. New developments and applications of the subsonic doublet-lattice method for nonplanar configurations//AGARD Symposium on unsteady aerodynamics for aeroelastic analysis of interfering surfaces, AGARD CP-80-71.— Nov. 1970.

6. Rodden W. P. The development of the doublet-lattice method/International forum on aeroelasticity and structural dynamics. — Rome. — 1997.

7. Мосунов B. A., Набиуллин Э. Н. Определение аэродинамических сил, действующих в дозвуковом потоке на упругие колеблющиеся поверхности, расположенные в разных плоскостях//Труды ЦАГИ. — 1981. Вып. 2118.

8. Appa K. Constant pressure panel method for supersonic unsteady airload analysis// J. of Aircraft. — Oct. 1987. Vol. 24.

9. Appa K., Smith M. J. C. Evaluation of the constant pressure panel method for supersonic unsteady airloads prediction//J.of Aircraft. — Nov. 1988. Vol. 26.

10. Кузьмин В. П., Ярошевский В. А. Оценка предельных отклонений фазовых координат динамической системы при случайных возмущениях.— М.: Наука.— 1995.

11. Евланов Л. Г., Константинов В. М. Системы со случайными параметрами. — М.: Наука.— 1976.

12. Kuzmin V. P. Estimation of wake vortex separation distances for approaching aircraft//Trudy TsAGI. — 1997. N 2627.

Рукопись поступила 21/XII2001 г.