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

Том XXXIII

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

№ 3—4

УДК 629.7.015.073

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

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

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

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

В литературе (см. [1—3]) принято рассматривать по отдельности недифференцируемые и дифференцируемые случайные процессы (в первом случае среднеквадратическое значение производной случайного процесса обращается в бесконечность). В общем случае вероятность превышения заданного большого порогового значения определяется формулой

P (х(т) > х* )«1 - exp (-t/T), (1)

где 0 < т < t, а значение T зависит от величины порогового значения и от характера процесса.

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

M[x(t )x(t + т)] = Kx (т) (2)

и спектральной плотностью

1

(®) = — [ С08 (®ТЖх (х^ X. (3)

2 % ■>

—да

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

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

j Sx(ю)ю

2d ю (4)

обращается в бесконечность, что характерно для моделей атмосферной турбулентности в форме Драйдена или Кармана. Для дифференцируемого случайного процесса задача облегчается благодаря возможности применения теоремы Крамера [1] к оценке среднего времени при больших значениях Я:

Т «л1- Кх^°-1 К (0)

ехр ( Я2Д) = л

ехр

( Я 72).

(5)

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

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

В простейшем случае для уравнения

х + х = 1^(0,

где Е>(^) — нормированный белый шум (тогда с х = 1), получим:

Т !

^ ехр (я2/2Уя .

(6)

(7)

В [3] приведены более точные значения параметра Т по сравнению с формулой (7) для нескольких значений Я, которые представлены в табл. 1 и могут быть использованы для контроля за точностью расчетов.

Т аблица 1

Я 2,5 3 3,5 4

точное 12,1 41,6 180 1007

приближенное 11,4 37,6 164 934

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

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

Т_1 =

(8)

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

0

Формулы, соответствующие более общему случаю недифференцируемого процесса, приведены в [2]. Если разность старших степеней ю знаменателя (m) и числителя (n) обозначить через m — n = 1 + а, где 0 < а < 2, то поведение безразмерной корреляционной функции процесса rx (т) = Кх (т)/Кх (0) при малых т определяется формулой

Гх (т) = 1 - C т |а + °^т|а), (9)

а формула, определяющая среднее время до первого пересечения гауссовским процессом одного из высоких уровней х = ± Rax, имеет вид

T *exp (rV2) R1-2/аС_1/^(а), (10)

где h(a) — функция от а, причем h(1) = 12, h(2) = ^%/2.

При а ^ 0 дисперсия процесса стремится к бесконечности. Значение C, согласно [3], определяется формулой

о % cos (%а/2)

C = -2 A2 —— -----------------> 0. (11)

sin (%(1 + а))Г(1 + а)

Значение A2 определяется из соотношения

A2 J Sx(ю)dю = A2 J Ют +... dю = 1, (12)

—ТО

—то

где

Sx (ю) = Sx (ю) lim

/" m—n \

ю

SXo)

(13)

Последний множитель в (13), очевидно, равен отношению коэффициентов при старших

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

Са в выражении (10) имеет размерность времени.

При а = 1 получим:

с = к*0)=—;—. (11.)

К (0) 7 _

1 1 | (ю) аю

—ТО

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

виде, отнесенном к времени т = —, где — масштаб турбулентности, и при допущении о

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

^г (®) =

л(1 + ш2 )

(14)

8мп (ю) = ■

1 + 3ю2

2я(1 + ю2 )

(15)

для продольной и поперечной составляющей соответственно. При этом для обеих составляющих справедливо соотношение

_2 ^

с мг, МП = 2

=1.

(16)

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

ч>п+ч>п+г =

г = <Щ(0, г + г = (/3-1^(Г)

(17)

где г — вспомогательная переменная.

Численные эксперименты проводились с разными шагами &, причем значение белого

шума на шаге задавалось постоянным, равным с/ ¡Л , где с — нормированная центрированная случайная величина с гауссовским распределением. Значения с на различных шагах некоррелированы друг с другом. Уравнения (6) и (17) интегрировались до тех пор, пока компоненты турбулентности не превышали по модулю последовательно значений Я = 2,5; 3,0; 3,5; 4. Начальные условия формировались случайным образом, в соответствии со

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

Р(Т) = ехр (-ТТ)/т .

(18)

Результаты расчетов значений Т, полученных при N = 200, Л£ = 0,001, приведены в табл. 2

для продольной (мг) и поперечной (мп ) составляющих.

Т аблица 2

Я 2,5 3 3,5 4

13,8 46,5 202 841

™п 9,4 32,2 120 660

Из сопоставления результатов, полученных в данном примере для компонента мг, и результатов, приведенных в табл. 1, следует, что численные расчеты при данном шаге и количестве реализаций не обеспечивают хорошей точности. Однако увеличить количество реализаций свыше определенного предела не удается из-за эффекта «зацикливания» датчика случайных (точ-

нее — псевдослучайных) чисел, особенно при Я > 4. В последнем случае требуемое коли-

2

0

чество независимых случайных чисел, распределенных по гауссовскому закону, составляет Іїї/At« 2 • 108, а требуемое количество независимых равномерно распределенных чисел превышает это значение примерно на один порядок. Анализ результатов показал, что последовательности случайных чисел в датчике при таком большом переборе начинают повторяться.

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

Аа

а„,

PTa

1 + PTa

(19)

не учитывающим изменения угла тангажа самолета, причем в качестве ориентировочного значения для безразмерной постоянной времени можно выбрать ^ = 0,4 (см. [6]). Поэтому система уравнений (17) была дополнена уравнением

й + а/Та = ™у,

(20)

описывающим изменение приращения угла атаки, и были проведены аналогичные численные эксперименты с фиксацией времени превышения этим приращением значений Raa, где ста = 0,623 (при ^ = 0,4), Я = 2,5; 3; 3,5; 4.

Результаты расчетов отражены в табл. 3.

Т аблица 3

Я 2,5 3 3,5 4

г 3,67 11,3 46,8 285

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

Т абл ица 4

Я 2,5 3 3,5 4

^а/ ^п 0,391 0,350 0,389 0,432

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

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

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

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

Ап

К

Аа 7аег р + Г

(21)

где Таег «—, Ь — средняя аэродинамическая хорда. Как указывалось выше, в этом случае

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

Т аблица 5

Т аег 0,01 0,02 0,04 0,06 0,1

°Апу 0,611 0,600 0,578 0,558 0,522

12,1 8,41 5,79 4,60 3,39

К Т (Я) Т (Я) Т (Я) Т (Я) Т (Я)

2,5 3,62 5,10 7,15 8,65 11,0

3 14,3 20,2 28,3 34.2 43,5

3,5 72,5 103 144 174 221

4 474 670 935 1090 1440

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

Таррг (Таег ) «^Ог + Т2 (Таег = °) • (22)

Здесь значения Теог определяются формулой (5) и содержатся в табл. 5, а значения Т (Таег = 0) определяются из табл. 3. Точность формулы (22) проверялась на основе сопоставления с численнымих результатами. Результаты сравнения отражены в табл. 6.

Т аблица 6

Формула (22) Таег = 0,01 0,02 0,04 0,06 0,1

К = 2,5 5,15 6,30 8,05 9,40 11,6

3 18,2 23,1 30,4 35,9 45,0

3,5 86,5 112 151 180 221

4 555 730 980 1120 1470

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

К = 2,5 5,45 6,35 7,95 9,85 11,3

3 17,1 21,5 29,0 32,8 45,0

3,5 81,0 107 158 190 229

4 565 785 1035 1120 1490

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

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

®шах «я/А• (23)

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

Рассмотрим для примера простейший недифференцируемый процесс вида (6). Потребуем, чтобы значение Т, определяемое по формуле (5)

T «ж

ж 2

я/Aí 2

* Ю ,

d ю

ехр

(я2Д), (24)

0 1 -ю2

оказалось в к раз меньше предельного значения (6).

В итоге получим:

А = к ГГ (25)

як 2 Я 2

(для к = 3 к < 0,004 при Я = 3, к < 0,002 при Я = 4).

В этом случае, в согласии с формулой (22), полученное в расчете среднее время окажется

завышенным всего лишь в Л +(г/к2 ) раз.

Как видно, с ростом К требования к допустимой величине шага ужесточаются. Рассматривая вышеприведенные примеры, можно убедиться, что в соответствии с формулой (10) отношения средних времен ТКУ!ТКХ и Та1 Тм>у составляют 2/3 и 1,5/3,4 = 0,44 соответственно, что подтверждается данными табл. 2 и 4 (для процесса изменения воздушного угла атаки отношение \К'х (0)|

-------- составляет 3,4).

Кх (0) , )

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

Эта спектральная плотность 8Х (ш) процесса х^) разделялась при ш> 0 на большое количество интервалов, вблизи центра каждого интервала выбиралась средняя частота шг-, случайный процесс аппроксимировался конечной суммой

(t) A sin (<V + Ф2 ). (26)

2=1

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

у^г +

А = 2 | ^(ш^о, (27)

“ _

где значения шг_, шг + являются границами интервалов, причем

“1_ = 0, ши+ = да.

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

с2 = 218Х (ш)а'ш . (28)

Далее вычислялись значения функции (26) с шагом Л£ = к до тех пор, пока при t = Т не выполнялось условие

\Ха ^^ > Ясх. (29)

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

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

ст Т

<30)

Интервал ±2с, соответствующий доверительной вероятности 0,955, при N = 400 составляет [0,9Т ^1,1Т ].

Очевидно, результаты численных экспериментов, проводимых для уровня «три сигма» (Я = 3), зависят также от величины шага к и от количества членов ряда п в формуле (26). Основная часть расчетов проводилась при использовании 48 и 54 членов ряда. Выбранные максимальные частоты шг для этих вариантов составляли соответственно 69,1 и 170. Совокупности частот формировались достаточно произвольно, но с учетом требования о том,

чтобы они не были пропорциональны членам натурального ряда: в противном случае процесс

(26) вырождается в периодическую функцию. Значение п должно быть достаточно велико, чтобы выполнялось неравенство:

П

£ А > тЯсх , (31)

г =1

где т > 2 ^ 3.

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

ТО

Е А2 = 2с2, (32)

г=1

эквивалентного соотношению (28), приводит к условию

0

n > -

m2 R 2

2

(33)

Таким образом, требуемое число членов ряда (26) возрастает с ростом Я. При этом шаг к должен удовлетворять условию

(34)

В конечном счете приходится выбирать шаг из условий компромисса между требованиями к точности расчетов и быстроте их проведения. Очевидно, что при возрастании порога (Я > 3,5 4) трудоемкость расчетов становится неприемлемой.

Спектральные плотности большей частью формировались таким образом, чтобы в

знаменателе содержались сомножители типа ^ш2 +1^ или ^ш2 +1^ , где I, т — небольшие целые числа.

По результатам расчетов можно считать, что при Я = 3, к < 0,005 точность расчетов является приемлемой. Об этом можно судить, сравнивая результаты расчетов для «эталонной» спектральной плотности (14) с точным решением (Т = 41,6 при Я = 3 ).

Однако использование столь малого шага заметно снижает скорость выполнения расчетов, тем более, что даже при к ^ 0 все равно сохраняется методическая ошибка, обусловленная конечным числом членов в формуле (26). Поэтому большое число расчетов было выполнено с шагом, существенно превышающим 0,005. Как показали результаты расчетов, при этом отношение значений параметра Т, вычисленных при различных шагах, мало изменяется при изменении спектральных плотностей. Тем самым с учетом значения Т = 41,6 для эталонной спектральной плотности (14) можно приближенно определять искомые значения Т при существенной экономии машинного времени. Зависимости вычисленных значений Т от шага к при п = 48 и 54 иллюстрируются в табл. 7.

Т аблица 7

h 0,05 0,02 0,005 0,005

n 48 48 48 54

N 800 400 1000 1000

T 59,34 50,5 48,34 44,63

Методическая погрешность, которая характеризуется коэффициентом завышения к(п, к) > 1, не исчезает при к ^ 0, причем, она убывает с ростом числа членов.

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

xa (0 *£[«« cos (<V) + v sin(<V)], (35)

2=1

<+

где щ, v — нормально распределенные числа с дисперсиями, пропорциональными J Sx (ю) dю .

<■Преимущество такой методики заключается в том, что случайный процесс является гауссовским при любом количестве членов ряда (35), в то время как процесс, описываемый формулой (26), приближается к гауссовскому при больших n согласно центральной предельной теореме.

Однако при использовании ряда (35) условие (28) и эквивалентное ему условие (32) выполняется неточно. Поэтому на больших интервалах времени (при больших R) при заданном пороговом уровне х* возникает разброс значения R = х*/ax, который уменьшается с ростом n

пропорционально 1/ П, но все же остается ощутимым даже при n > 25 + 50.

П

Учитывая, что функции, определяющие зависимость от значения R среднего времени, T(R) * C exp(Д2/2) или T(R) * C exp(r2/2)/R, имеют большую вторую производную T'(R),

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

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

Sx (ш) =

ш

(ш2 + a2)(ш2 + (2 - а)2)

(36)

В данном случае, используя формулу (11'), убедимся, что значение С не зависит от значения а. Результаты расчетов при Я = 3, п = 54, к = 0,005, N = 200 приведены в табл. 8.

Т аблица 8

а 0 0,4 0,5 0,7 0,8 0,9 1,0

T 25,73 26,42 24,70 26,43 26,54 23,31 25,20

Учитывая формулу (30) для разброса результатов расчета, можно сделать вывод о том, что для спектральной плотности вида (36) соотношение (11') подтверждается.

В разделе 1 рассматривался вариант, когда дифференцируемый случайный процесс при Тег ^ 0 становился «почти недифференцируемым», и для этого варианта была предложена полуэмпирическая формула (22).

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

S (ш) = aSd (ш) + bSn (ш),

(38)

где Sd (ш), Sn (ш) — спектральные плотности дифференцируемого и недифференцируемого процесса; a, b — коэффициенты; значение R не очень велико. Отметим, что такое представление достигается не единственным образом. Сложность поставленной задачи усугубляется тем, что при больших R для недифференцируемого (при а = 1) и дифференцируемого процессов среднее

время пропорционально exp (R2l2)jR и exp (R2/2) соответственно.

Представим простую аппроксимационную формулу в виде дробно-рациональной функции от коэффициентов а и b, считая коэффициент к(n, h) неизменным:

T * к(h, n) -

—a + -0b

а + -3bR

exp ( R2 /2).

(39)

При Ь ^ 0 формула (39) определяет значение Т для дифференцируемого процесса. Отсюда, пользуясь теоремой Крамера [1], получим

же

Sxd (ra)d ш

- = -

xd

xd

j<jy2Sxd (со)б/со

о

При больших значениях R, используя соотношение (11'), найдем:

\ = Я Kxd(°) ^2 = RKn

3 b |KXn(0)Г \K'„

(0)_

(0)1'

Применимость предложенной формулы проверялась при п = 48, к = 0,02, N = 1000 для случая

(ш) =

(l + ш2 ) 1 ■

ш

2 '

(42)

Коэффициент завышения определялся путем сопоставления результатов расчета для спектральной плотности lj(l + ш2) (T = 308) с теоретическим результатом (применение

теоремы Крамера, T * % exp (R2Д) = 282) . Формула (39) в этом случае приобретает вид:

a + 2b

T * 308-

a + 15b

(43)

Если отказаться от предположения о неизменном коэффициенте завышения и добиваться совмещения результатов расчета в «реперных» точках (а = 1, Ь = 0 и а = 0, Ь = 1, последней точке соответствует по результатам расчета значение Т = 50,5), то значение е3 следует изменить, и модифицированная формула (43) приобретает вид:

T * 308-

a + 2b a +12,2b

(44)

Сравнение результатов численного расчета с формулами (43) и (44) дано в табл. 9.

Т аблица 9

b (a =1) T Формула (43) Формула (44)

0 308 308 308

0,05 191,2 193,6 210,4

0,1 151,9 147,8 166,5

0,2 118,5 107,8 125,3

0,3 102,1 89,6 105,8

0,5 80,2 72,5 86,8

1 66,3 55,3 70,0

2 59,3 49,7 60,6

3 56,1 46,9 57,3

5 53,7 44,6 54,6

12 52,5 42,5 52,2

20 51,5 42 51,5

a = 0 (b =1) 50,5 41,1 50,5

a = -1/3 44,4 35,1 43,3

a = -2/3 37,5 28,6 35,5

a —1 29,8 22,0 27,5

Аналогичные расчеты были проведены для спектральной плотности вида

а Ь

Ях (®) =

(1 + ш2)3 1 +

2 '

(45)

ш

В предположении к(n, h) = const формула (39) приобретает вид:

.a + 2,67b

T * 545

a + 34,6b ’

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

T * 545

a +2,67b

а + 28,8й

Сравнение результатов расчета с выписанными формулами при Я = 3 дано в табл. 10.

Таблица 10

(47)

a b T Формула (46) Формула (47)

1 0 545 545 545

1 0,05 242 226 253

1 0,1 170,6 154,8 178

1 0,2 117,9 105,5 123,6

Продолжение табл. 10

a b T Формула (46) Формула (47)

1 0,3 99,6 86,2 101,8

1 0,5 77,6 69,5 82,6

1 1 65,3 49,2 67,0

1 2 58,2 49,2 58,9

1 3 54,9 46,8 56,1

1 5 53,6 44,9 53,9

1 12 52,5 43,2 51,9

1 20 51,4 42,7 51,3

0 1 50,5 42,0 50,5

-1/3 1 45,3 37,2 44,7

-2/3 1 40,5 32,2 38,8

-1 1 34,5 27,1 32,7

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

Модифицированную формулу типа (44), (47) можно применить и к случаю, когда спектральные плотности не описываются дробно-рациональными функциями частоты. Пусть спектральная плотность имеет вид:

(®) = 7

■ + ■

~-|11/6 Г 9“|5/6'

1 + (1,339ш)2 1 + (1,339ш)2

(48)

Возьмем за основу результаты расчетов при п = 35, И = 0,02, Т = 528 при а = 1, Ь = 0; Т = 35,3 при а = 0, Ь = 1. Тогда формула, подобная (44), записывается в виде:

T * 528

a + 2,5b

а +37,4Ь

Результаты расчета сравниваются с этой формулой в табл. 11.

a b N T Формула (49)

1 0 300 528 528

11 1 200 157 147

5 1 400 90 93,4

2 1 400 60 60,3

0 1 800 35,3 35,3*

-5 8 800 27 26,9**

-1 1 800 24,1 21,8

*, ** — спектральные плотности продольных и

поперечных порывов для модели Кармана.

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

ЛИТЕРАТУРА

1. Стратонович Р. Л. Избранные вопросы теории флюктуаций в радиотехнике. —

М.: Советское радио.— 1961.

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

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

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

5. Swolinsky M. Turbulence engineering models, aircraft response//Report of Technical University Braunschweig.— 1995.

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

№ 1 — 2.

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

Рукопись поступила 19/II2001 г.