Применение вероятностных методов к задачам динамики полета Текст научной статьи по специальности «Математика»

Том XXXIX

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

№ 1 — 2

УДК 629.735.33.015

ПРИМЕНЕНИЕ ВЕРОЯТНОСТНЫХ МЕТОДОВ К ЗАДАЧАМ ДИНАМИКИ ПОЛЕТА

А. В. БОБЫЛЕВ, В. А. ЯРОШЕВСКИЙ

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

1. За последние годы задача обеспечения безопасности полетов приобретает все большую актуальность. Это выражается в формулировании требований об обеспечении благополучного исхода полетов с высокой степенью надежности. То есть, допустимые вероятности возникновения тяжелых летных происшествий выражаются числами 10-6 —10-10 на час или один полет. В связи с этим появляется необходимость моделирования процесса полета в очень большом объеме с учетом всех существенных случайных возмущений.

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

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

x(t)s Z Ck cos((Bkt) + dk sin((Bkt), (1)

k=1

где Ck, dk — случайные коэффициенты.

Обсудим те трудности, которые возникают при моделировании случайного процесса в виде суммы гармоник. Методика, описанная в [1], заключается в том, что заданная спектральная плотность Sx (ю) процесса x(t) разделяется при ю>0 на большое количество интервалов, вблизи

центра каждого интервала выбирается средняя частота юг-, и случайный процесс аппроксимируется конечной суммой

П

(2)

Коэффициенты 4 определяются по формуле:

(3)

где значения юг-_, юг+ являются границами интервалов, причем ю1- = 0 , юп+ = ^ . Фазы фг задаются с помощью датчика случайных чисел, равномерно распределенных в интервале [0, 2п]. Дисперсия процесса определяется по формуле:

Процесс (2) при большом числе членов ряда близок к гауссовскому согласно центральной предельной теореме. Далее вычисляются значения функции (2) с шагом А( = к до тех пор, пока при ^ = Т не выполняется условие

где Я — «большое» число. На практике это число считается большим, если оно не меньше 3. Преимущества представления случайного процесса в виде (2) по сравнению с (1), где ек, ёк — случайные коэффициенты с гауссовским распределением, обсуждаются в работе [1].

При достижении условия (5) значение Т запоминается, задается новая случайная совокупность фаз ф1 фп, и расчет повторяется, пока число реализаций не достигает заданного значения N.

Параметр Т определяется путем усреднения значений Т. Распределение значений Т является экспоненциальным, при этом среднеквадратическое отклонение оценки параметра Т в результате выполнения N экспериментов составляет т/тЩ . Таким образом, например, при числе экспериментов, равном 1000, среднеквадратическая погрешность определения среднего времени составляет около 3%. Поскольку описанный процесс вычисления среднего времени является очень трудоемким и трудоемкость сильно возрастает с ростом Я, ограничимся в данной статье представлением численных результатов только для случая Я = 3. Последовательность частот задавалась произвольно. Максимальные значения частот зависели от числа членов ряда (2): при п = 54, 62, 67, 75

эти значения юп_ составляли 153, 273, 358 и 478 соответственно. Отметим, что последовательность частот не должна быть пропорциональна членам натурального ряда, во избежание реализации периодического процесса. С ростом Я потребное число членов суммы (2) возрастает.

2. Рассмотрим простейший пример применения указанной методики. При анализе динамики полета в условиях атмосферной турбулентности часто используют формулы Кармана, описывающие изменение по времени атмосферной турбулентности для продольных и поперечных порывов, [2]:

(4)

0

\ха ()|^ Я°х,

(5)

1

(6)

* (ю) = 1 + (8/3)(1-34ю)

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

- У

ному времени, равному ^ =—, где — масштаб турбулентности. Результаты расчетов среднего времени Тдля случая Я = 3, N = 1000 при разном числе слагаемых в сумме (2) и при разных значениях шага (в том числе — при большом значении, вторая строка, [1]) приведены в первом и втором столбцах табл. 1. Тем самым по результатам, приведенным во втором столбце, можно оценить частоту выбросов значения ветрового угла атаки а№. Результаты, приведенные в треть-

2// 2\п/6

столбце, относятся к случаю 51(ю) = ю/ (1 + (1.34ю) ) .

Таблица 1

ем

Б (ю) т = 0.2 т = 0.4 т = 0.6 п Н

35.3 27.0 24.1 35 0.02

25.3 18.1 5.93 7.84 10.25 54 0.005

21.5 15.2 10.8 3.54 5.71 6.65 75 0.002

20.3 15.1 11.5 75 0.001

Видно, что в данном случае уменьшение шага ниже значения 0.002 нецелесообразно, учитывая, что трудоемкость расчетов пропорциональна TnN|h.

Вместе с тем, специалистов в области динамики полета чаще всего интересуют возможные «выбросы» значений поперечной перегрузки, которые примерно пропорциональны суммарному углу атаки а + а№ (или нелинейной функции от этого угла). Поэтому следует учесть, что указанные углы связаны передаточной функцией, которую (для безразмерного времени) можно в пер-

Тр

вом приближении аппроксимировать выражением

Тр +1

где р — оператор дифференцирования,

т = 0.2 -0.6, [3]. В трех последующих столбцах приведены значения среднего времени, соответствующие спектральной плотности Бп, умноженной на квадрат модуля передаточной функции

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

тУ 1 + т2ю2

выбросов ветрового угла атаки. Этот факт следует иметь в виду при определении в реальном случае реакции самолета на расчетный порыв ветра. Действительно, в некоторых случаях принято определять величину порыва ветра, соответствующую заданной вероятности ее превышения (определяемой значением Я) за заданное время полета [4, 5], затем анализировать реакцию самолета на этот порыв и вычислять максимальную величину полного угла атаки. На самом деле, предельный суммарный угол атаки, полученный в результате этого расчета, может достигаться в 2 — 3 раза чаще.

3. Как видно из предыдущего примера, спектральная плотность процесса не является дробно-рациональной функцией частоты. Этому случаю посвящены, в частности работы [6, 7]. Пусть спектральная плотность процесса определяется формулой:

, -/7 ,

ю +... і

ю +...

= Л2 Б (ю),

(8)

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

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

К (т) = 1 - ("| т|“ +...,

(9)

где

с =

Кх

(0),

а = т - п -1.

Зададим простейшее семейство спектральных плотностей в виде:

Б (ю) =

ю

21

(1 + ю2)+1

где к = {т -п)/2 > 1/2 (при к < 1/2 дисперсия процесса обращается в бесконечность). Для того

чтобы проиллюстрировать зависимость результатов численного расчета от шага и количества членов в сумме (2), приведем некоторые из этих результатов для случая I = 0, Я = 3 в табл. 2

Т аблица 2

к п Н N Т к п Н N Т

0.7 54 0.005 400 10.6 1 54 0.002 900 44.8

0.7 54 0.002 2000 10.5 1 57 0.005 500 44.7

0.7 62 0.005 800 9.05 1 62 0.005 6400 41.6

0.7 62 0.002 1000 8.68 1 62 0.002 1600 41.3

0.7 67 0.002 1600 7.55 1 75 0.005 1000 41.8

0.7 67 0.001 1600 7.25 1.5 54 0.002 200 170

0.7 75 0.001 1000 6.66 1.5 54 0.005 800 174

1 54 0.05 1000 63.2 1.5 62 0.005 1200 163

1 54 0.02 800 51.6 2 54 0.01 1200 315

1 54 0.005 700 45.1 2 62 0.01 1350 329

Среднее время, полученное в численных экспериментах, почти всегда оказывается завышенным по сравнению с истинным. Это время уменьшается с увеличением числа членов ряда и с уменьшением шага. Естественно, что наиболее жесткие требования к этим параметрам предъявляются при малых значениях к, поскольку при больших к спектральная плотность убывает быстрее по мере роста частоты. В дальнейшем будем считать наиболее достоверными результаты, соответствующие максимальным значениям п и минимальным значениям h. Отметим, что для классического случая к = 1 существует точное решение Т = 41.6, [2], и результаты численного эксперимента при п > 62 оказываются близкими к теоретическому значению (впрочем, это может объясняться простым совпадением).

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

Т = і (Я, к)ехр(72).

(11)

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

конечно, а для недифференцируемого ( > 3/2) — бесконечно. Основная трудность заключается в том, чтобы среди всех пересечений выделить первые пересечения уровня х = ±Ясх . При больших значениях Я допустимо отождествить первые пересечения со всеми пересечениями. Тогда функция t для недифференцируемого процесса определяется формулой [7]:

1

- Я2 (I - Кх (і))■ ^ткхй I 2(+Кх())

Здесь Кх ^) — корреляционная функция нормированного случайного процесса (х (0) = 1). При очень больших значениях Я из этой формулы получим асимптотическое выражение

> у/(2к-1)

а = с (к) | Б (ю)Ою/Я2

V 0

где функция С (к) определяется формулой:

Я,

(13)

5-2 к

С (к) = 22(2к-1)г

3 - 2к V Г(к)( -1/2)^(2 к-1)

2(2к -1) ^^л/Пг(3/2 - к) Значения С (к) приведены в таблице 3.

'(22 -1).

(14)

Таблица 3

к 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4

с (к) 2.02 0.898 0.756 0.749 0.798 0.894 1.05 1.34 2

Формула (13) удобна тем, что значения коэффициента С (к) имеют порядок единицы и слабо изменяются в широком диапазоне значений к. Для случая I = 0 из этой формулы получим результат, совпадающий с результатом, приведенным в работе [7]:

3-2к

іа:, =(1/Я)2к-1 22(2к-1) Г

3 - 2к " Г( +1/2 )"

2(2к -1) Г (3/2 - к)

V (2к-1)

/(2к -1) при к < 3/2.

(15)

Формулу (13) можно объяснить из соображений о том, что среднее время должно быть пропорционально характерному времени корреляции, которое, с учетом формулы (9), пропорцио-

/1/ г \1(2к-1) нально 11/ С I .

Для дифференцируемого процесса известная асимптотическая формула приобретает простой вид:

|Б(ю)Ою/ | ю2Б(ю)Ою при к >3/2.

(16)

В случае (10) при I = 0 получим простую формулу = П2( -3/2). Можно убедиться,

что при заданном значении Я функция tas по мере увеличения к от 0.8 до 1.5 возрастает и стремится к бесконечности пропорционально 1Д/3/2 — к, а затем обращается в нуль при к = 1.5 и возрастает при увеличении к. В то же время очевидно, что в реальности в окрестности к = 1.5 эта функция должна изменяться непрерывно. Поэтому при выполнении численных экспериментов этой окрестности было уделено особое внимание. В табл. 4 результаты численного расчета среднего времени сравниваются с результатами расчетов по приближенным формулам (12), (13) и (16).

Таблица 4

к 1 Численный расчет Формулы 13, 16 Формула 12 п Н N

0.7 0 6.66 17.4 9.07 75 0.001 1000

0.7 1 3.92 6.20 2.89 75 0.001 1000

0.7 2 3.26 4.13 2.12 75 0.001 1000

0.7 3 2.80 3.19 1.76 75 0.001 1000

к 1 Численный расчет Формулы 13, 16 Формула 12 п h N

0.8 0 13.7 20.7 17.6 75 0.002 1000

0.8 1 7.12 8.83 6.44 75 0.002 1000

0.8 2 6.07 6.28 4.76 75 0.002 1000

0.8 3 4.87 5.05 3.95 75 0.002 1000

5/6 0 16.5 75 0.002 1000

5/6 1 8.93 75 0.002 1000

0.9 0 25.5 27.8 28.7 75 0.002 1000

0.9 1 13.6 13.0 11.5 75 0.002 1000

0.9 2 10.3 9.51 8.55 75 0.002 1000

0.9 3 8.81 7.80 7.10 75 0.002 1000

1.0 0 41.6 37.6 42.9 62 0.005 6400

1.0 1 22.9 18.6 18.4 62 0.005 1000

1.0 2 18.7 13.8 13.8 62 0.005 1000

1.0 3 16.3 11.5 11.5 62 0.005 1000

1.1 0 61.6 50.8 61.8 62 0.005 800

1.25 0 104 82.6 107 62 0.005 1150

1.25 1 52.9 44.8 51.6 62 0.005 800

1.25 2 38.7 34.1 39.1 62 0.005 800

1.25 3 35.5 28.6 32.7 62 0.005 800

1.4 0 139 165 190 62 0.005 800

1.45 0 161 255 62 0.005 1000

1.5 0 163 62 0.005 1200

1.5 1 88.5 62 0.005 600

1.5 2 70.8 62 0.005 600

1.5 3 60.1 62 0.005 600

1.6 0 198 127 62 0.005 1200

1.75 0 240 200 62 0.005 800

2.0 0 329 283 62 0.01 1350

2.0 1 174 163 62 0.01 1000

2.0 2 126 126 62 0.01 1000

2.0 3 121 107 62 0.01 1000

Как видно, при Я = 3 значения, вычисленные по формулам (12), (13), близки к результатам расчета при 0.9 < к < 1.25 .

Изобразим на рис. 1 — 4 полученные результаты в виде точек, определяющих безразмерное время t, а также кривые (12) (жирные линии), (13) и (16) в функции от параметра к при Я =3, I = 0, 1, 2, 3. Напрашивается вывод о том, что асимптотические кривые (13) и (16) проще всего «срастить», проводя к ним общую касательную. Этот вывод также можно сделать, исходя из результатов, полученных В. П. Кузьминым в работе [7], основанных на статистическом моделировании в большом объеме по специальной методике.

Рассмотрим случай I = 0 и определим значения t с помощью графических построений. Убедимся в том, что при использовании этой гипотезы в окрестности к = 3/2 значения t приближаются к асимптотическим (формулы (13), (16)) очень медленно с ростом Я. В табл. 5 приведены результаты вычислений.

Таблица 5

Я 5 20 50 200 1000 5000 25 000

7 (3/2) 1.56 1.4 1.32 1.22 1.13 1.05 1

Можно показать, что при очень больших Я среднее время примерно пропорционально 1/л/1п Я. На рис. 5 приведена зависимость значения ? (3/2) от 1п Я. Как и следовало ожидать, в пограничной точке, при к = 3/2, зависимость является логарифмической.

/

/ /

-2 5- • /

? / у

■ • • /

1. Э / /

" 1 " . Л ч. /

и. 3 л У г

0,8 1 1.2 1,4 1,6 1,8 2 к

Рис. 1. Бх (<в) = 1 (1 + ю2 ) :

— формула (12);

— формулы (13), (16);

— численный расчет

1 9^-

1 .

1 -П 75 - 1 • '

и. / 3 Р1 5 -

и. 3 / /

"1).2У У

0.8 1 1.2 1.4 1.6 1.8 2 к

Рис. 3. Бх(ю) = о>У(1 + ю2)+2 :

— формула (12);

— формулы (13), (16);

— численный расчет

/

.... 1 *ч- / V

1 _

П 5 •

и. э У У

0.8 1 1.2 1,4 1.6 1.8 2 к

Рис. 2. Бх (ю) = о>У ( + ю2) :

— формула (12);

— формулы (13), (16);

— численный расчет

1 1 2

1 /

/ 1

1).о 0 6 ! • /

0 4 ! '

0 2 у 1/

/•

0.8 1 1.2 1.4 1.6 1.8 2 к

Рис. 4. Бх (ю) = а>У ( + ю2) :

— формула (12);

— формулы (13), (16);

— численный расчет

Рис. 5. k = 3/2

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

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

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

плотности случайного процесса. Отметим, что для спектральной плотности вида А&Рj(2b2 +1)

по сравнению с А&Рj(2 +1) среднее время просто умножается на b, что согласуется с асимптотической формулой (13). Например, сравнивая результаты расчета для случая радиальных порывов при k = 5/6, приведенные в табл. 1 и 4, получим, что среднее время в табл. 1 превосходит результат в табл. 4 в 1.3 раза (по сравнению с 1.34).

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

Sx (ю) =

А1ю'

А2Ю”2

+... +

A a"k

(t2+1) (2t22+1)2 (m+1)

(17)

Значения п, в частности, могут быть равны нулю. Выделим наименьшее из значений 2кг- — пI. Оно определяет асимптотическое выражение для среднего времени при достаточно больших Я . Рассмотрим простейший случай, когда сумма (17) включает два слагаемых и значение 2кг — п является наименьшим. Тогда для вычисления среднего времени можно применить асимптотическую формулу (13), в которую следует подставить выражения:

k=k-V2, s(ю)=SxИт^Д.

(18)

Но если значение А2/?22к2 заметно превосходит А^Т12к , то при не слишком больших Я среднее время будет определяться в основном вторым слагаемым и может сильно отличаться от асимптотического значения (что особенно ярко проявляется, например, при к ~3/2, I = 0). Тогда можно предложить полуэмпирический метод уточнения среднего времени, соответствующего сумме вида (17), состоящей из двух слагаемых, если для каждого из этих слагаемых среднее время, соответствующее данному не слишком большому значению Я, известно. Если предположить, что асимптотическая формула (13) применима, то безразмерное среднее время определяется формулой:

J Sx (<o)d _0_________

(T12k1 W

Таким образом, числитель включает слагаемые, пропорциональные как Ау, так и А2, а знаменатель не зависит от А2.

Если эта формула при умеренных значениях Я и больших А2 недостаточно точна, то можно включить в знаменатель слагаемое, пропорциональное А2 таким образом, чтобы при Ау = 0 и А2 = 0 значение среднего времени совпало с результатом численного расчета. В итоге получим эмпирическую формулу:

Здесь спектральные плотности Бху(ю), Бх2(ю) соответствуют первому и второму слагаемому в сумме вида (17), значения Туп, 72„ определяют среднее время, полученное на основании численных экспериментов при заданном значении Я. Подобный метод применялся и в [1] для случая к = 1. Отметим, что эта формула согласуется с асимптотической формулой (13) при 2ку — п = 2к2 — п2. Проверим эту методику на примерах:

(по сравнению с 6.66);

в) определим средние времена для модели Кармана, исходя из данных, приведенных в табл. 1. В данном случае к = 5/6. Для радиальных компонентов турбулентности имеем Тг = 21.5. Для

зуем формулу (20) и получим Т - 16.3 (по сравнению с 15.2).

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

1. Бобылев А. В., Ярошевский В. А. Методика проведения численных экспериментов для решения задач о больших выбросах случайных процессов // Ученые записки ЦАГИ. 2002. Т. XXXIII, № 3 — 4.

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

3. Парышева Г. В., Ярошевский В. А. Проблема формирования расчетных ветровых возмущений для задач динамики полета // Ученые записки ЦАГИ. 2001. Т. ХХХ11, № 1 — 2.

4. Проект отраслевого стандарта «Модель турбулентности атмосферы». 1984.

5. Селихов А. Ф., Чижов В. М. Вероятностные методы в расчетах прочности самолета. — М.: Машиностроение, 1987.

6. Лидбеттер М., Линдгрен Г., Рутсен Х. Экстремумы стационарных случайных последовательностей и процессов. — М.: Мир, 1989.

7. Кузьмин В. П. Метод оценки вероятностей больших выбросов случайных процессов // Ученые записки ЦАГИ. 2003. Т. ХхХ1У, № 3 — 4.

(

Л1/ (2к—1)

| Бх1(ю^ ю+15х 2 (ю)а? ю

Т -

0

0

(20)

|Бх1 (юУ®/Т12пк—1 + |^х2 (ЮУ®/Т2

т2 к—1

2п

0

/

Т1п = 329, Т2п = 22.9, Т - 42.8 (по сравне-

нию с 41.6);

к = 0.7, Т1п = 225, Т2п = 3.92, Т - 7.52

ЛИТЕРАТУРА

Рукопись поступила 7/112007 г.