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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том VII 1976

№ 4

УДК 629.735.3.018.4:620.178.3

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

Г. В. Вронский

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

Определение усталостного повреждения в расчетных точках авиационной конструкции является центральной задачей расчета ресурса планера на этапе проектирования самолета. В настоящее время общее признание получил метод расчета повреждаемости по полным циклам [1]. Известно, что с помощью метода Монте-Карло [2] можно получить ансамбль возможных временных реализаций, соответствующих заданной спектральной плотности мощности стационарного эргодического случайного воздействия со средним значением, равным нулю.

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

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

Для преодоления этих трудностей предложено несколько гипотез, с помощью которых при дополнительных предположениях относительно закона распределения нагрузок можно было получить оценки ресурса сверху и снизу [4, 5]. В. Л. Райхером [6] из некоторых энергетических соображений была предложена „гипотеза спектрального суммирования“, позволяющая свести нагружение с непрерывным спектром к некоторому эквивалентному в смысле повреждаемости гармоническому воздействию. При этом статистика полных циклов не используется для расчета повреждаемости.

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

Предлагаемый алгоритм запрограммирован на машинном языке БЕМШ для ЭЦВМ БЭСМ-6 и оформлен в виде модуля комплексного прочностного расчета. Данная программа достаточно универсальна, так как оперирует с заданными амплитудно-фазовыми частотными характеристиками и обладает эффективной скоростью счета, позволяющей формировать достаточно представительные статистики. В частности, для ансамбля из 100 входных реализаций длиной в 512 отсчетов по времени расчет повторяемости нагрузок и повреждаемости, включающий учет всех шести компонентов напряжений, занимает 3—4 мин. В качестве примера использования данного метода приведен расчет повреждаемости элемента конструкции при полете летательного аппарата в неспокойной атмосфере.

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

о ! \ 2 ^ 1 З£3 /1 ч

о„, (со) == о2-------. (IV

•*£/ (1 -Н2)2

Здесь — среднеквадратичная скорость порыва;

Ь — масштаб турбулентности; и—скорость полета; со — круговая частота; безразмерный параметр £ выражается формулой

__о

= И'

Количество вычисляемых спектральных ЛИНИЙ Sw(w) определяется числом заданных спектральных линий частотной характеристики. Дисперсия Dw вертикальных порывов соответственно будет:

со

Dw = j *£«,(“) da.

о

Будем считать, что частотные характеристики, соответствующие расчетной точке самолетной конструкции Ra/(iсо), заданы таблично. Таким образом, задана комплексная матрица [/?(го>)], имеющая шесть строк, в которых с некоторым заданным шагом Дсо по круговой частоте со записаны частотные характеристики для напряжений ах, <зу, <зг, хху, xxz, xyz соответственно. Число столбцов этой матрицы, которое определяется числом отсчетов по со значений частотных характеристик, обозначим через Л^.

Обычно полагают, что скорость вертикальных порывов W (t) является стационарным нормальным эргодическим случайным процессом [8]. Тогда на конечном интервале времени [О, Т\ можно написать ее каноническое разложение [2]:

со

W (t) = 2 SR (k) cos (kA соt) -f St (ft) sin (йД соt). (2)

*=i

Здесь период T определяется шагом по частоте Дсо, с которым дискретизируется спектральная плотность дисперсий порывов:

Дш

Согласно методу канонических разложений [2], 5«(й) и 5/(&) представляют соответственно действительную и мнимую части спектра амплитуд вертикальных порывов и являются некоррелированными случайными величинами с математическими ожиданиями, равными нулю. Поскольку ІІ7(¿) является нормально распределенной случайной величиной, £#(&) и 5/(6) будут также распределены нормально [2].

Дисперсии величин .$#(&) и 5/(А) одинаковы для каждой пары с одним и тем же индексом & и равны дисперсии порывов, приходящейся на интервал частот Дсо в окрестности со, т. е. равны (со) Дю:

Е> [ЗД)1 = О [5/ {к)] = (¿Д«>) Дсо.

Очевидно, что случайная функция \№(і) обладает нулевым

средним значением. Используя формулы Эйлера для комплексного

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

СО

£ ф»*е*км>

й = —со

где і — мнимая единица, а комплексные коэффициенты Фурье Ф^, выражаемые формулой

Л _ 5я(й) — (ъ\

— —~ ^ ^

представляют собой амплитудный спектр вертикальных порывов квазидетерминированной периодической функции №{$) на выбранном периоде разложения Т.

Итак, по спектральной плотности дисперсий порывов, выражаемой формулой (1), может быть получен ансамбль из N квазидетер-минированных реализаций №(() на интервалах времени длиной Т. В каждой из этих реализаций случайные величины 5^((о) и 5/(<и) выбираются методом статистических испытаний (Монте-Карло) из датчика нормальных случайных чисел. Ансамбль N реализаций (./V—достаточно большое целое число) скорости вертикальных порывов на интервалах времени Т может считаться приближением к стационарному случайному процессу с непрерывным спектром дисперсий (1).

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

Здесь — частотная характеристика /-го компонента на-

пряжений, определенная для положительных частот и задаваемая матрицей Фщ,(гш) — спектр комплексных амплитуд вертикаль-

ных порывов, определяемый формулой (3).

С помощью формулы (4) определяем комплексный спектр напряжений при положительных частотах. Применяя быстрое преобразование Фурье к спектральному разложению компонентов напряжений, и используя свойство вещественности компонентов напряжений [9], получим на отрезке времени [О, Т] временные реализации каждого компонента динамических напряжений в выбранной точке конструкции, соответствующие данной реализации И^(0:

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

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

Определим вначале главные напряжения из характеристического уравнения тензора напряжений [10]:

где значения инвариантов напряженного состояния /г (г=1, 2, 3) задаются формулами

= /?<гу(го>) Фю (ш).

(4)

<з3 — ^ О2 — /2 а — У3 = 0,

(5)

'ху *хг "угщ

Корни уравнения (5) -г- главные напряжения расположим в порядке невозрастания величин:

а1 ^ 32 аз-

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

а, = ( 01 ПрИ 1а11>1аз1'-

1 1 а3 при | о3 | > | б!

Введем еще удобное для дальнейшего изложения обозначение:

о3 при \в11 > I о31;

аш ,

а1 при | 38 I > I 3!

С учетом знака напряжений [3] в качестве эквивалентных напряжений в смысле усталостного повреждения будем считать:

I теория (максимальные по модулю главные напряжения):

аэкв °1>

II теория (максимальные по модулю деформации):

°экв = °i — И- (а2 + аш);

III теория (максимальные по модулю касательные напряжения):

^ЭКВ --- ®IIli

IV теория (интенсивность напряжений по Мизесу, которой припишем знак изменения объема элемента ДУ):

39KB = °HSign(AI/). (6)

Здесь ои — интенсивность нормальных напряжений, получающаяся из выражения энергии формоизменения [10]:

3И — У(^ — оу)2 + (оу — зг)2 + (в, - °х)2 + б (*ху + х2г + . (7)

Изменение объема можно выразить через напряжения [11]:

AV=(ex + zy+ ег) у=(ох + ву + 0e)Z(l _2[а), (8)

Г.

где гх, eyt ez — относительные деформации; V — объем элемента; Е— модуль Юнга; ^ — коэффициент Пуассона.

Учитывая, что значение jx не превышает 0,5, из выражений (6) — (8) окончательно получим:

°IV экв = °и Sign (рх — Оу + <32) = Ои Sign Л.

Вычислив временную реализацию эквивалентных напряжений o3KB(/f) по одной из четырех указанных теорий (выбор теории зависит от материала и вида напряженного состояния, для большинства пластичных материалов в лучшем соответствии с экспериментом находится теория энергии формоизмерения [12]), будем использовать ее при подсчете повреждаемости за время Т. Выделив из нее экстремумы и обработав их с помощью метода полных циклов [1]. получим таблицу повторяемости полных циклов за время Т с двумя параметрами — средним значением напряжений цикла ет и

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

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

а-1 р„.

к=0 /=о 14 рк!

Здесь Мт — число интервалов, на которые разбит диапазон изменения а,„; Ыа — число интервалов, на которые разбит диапазон изменения оа\ Рк] — повторяемость полных циклов, попадающих в А-й диапазон по оя и у'-й диапазон по <зт:

¿ = 0, 1, . . ., (№а — 1);

7 = 0, 1, . . ., (Л^-1);

¿А'за<оа<(А+ 1)

(ат)т1п “I- ^ (°т)тт Ч~ (/ ~Ь 1) А<3т’

где (от)т{П — минимальное значение ат; Дзот — шаг по от; Аоа-—шаг по оа.

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

2. В качестве примера приведем результаты расчета напряжений, повторяемости полных циклов и повреждаемости в расчетной точке элемента крыла с концентратором, амплитудно-частотная характеристика |/?(го>)| которого, заданная в 179 точках, приведена на фиг. 2. Считалось, что обшивка изготовлена из сплава Д1, при этом использовалась кривая выносливости с запасом по напряжениям Ч» = 1,3.

Расчет проводился при нулевой статической нагрузке и следующих значениях параметров: ¿7= 222 м/с; Д<о = 0,5 рад/с; от варьировалось от 5 до 15 м/с, Ь — от 20 до 300 м.

Число N разыгрываемых методом Монте-Карло случайных временных реализаций 47 (¿) равнялось 100.

В табл. 1 приводится повторяемость полных циклов за время Т= 12,56 с при ода==5м/с; ¿ = 20 м.

В табл. 2 приведены значения повреждаемости П за время Г= 12,56 с как функции параметров Ь и aw. Как показали поверочные расчеты, при фиксации одного из этих параметров зависимость 1дП от логарифма другого параметра приблизительно линейна, в

§

Фиг. 1

10

0,5

\<à)\

,

1

AE<3æ

1

I \«х

V У V J 1

і \ 1 ^ ¥ — \ \ \ V

J 1 \ /' J

/ '!г V І=Г — -— ;С S • \

0 10 20 Л7 Í/7 50 60 70 ЗО о>} рад/с Фиг. 2

L, м 20 ЗО 40 43 ^

5 0,367-10-8 0,836-10-1° 0,206.10-ю 0,12Ы0-Ю

7 0,129-10“7 0,148-10-* 0,469-10-9 0,294-10-8

8 0,472-10-7 0,606-10-8 0,239-10-8 0,874-10-9

10 0,278-10-6 0,785-10—7 0,31-10-7 0,164-10-7

13 0,168-10-5 0,341-10-е 0,123-10-6 0,11-10-6

15 0,41-10-5 0,148-10-5 0,566- 10-е 0,487-10-6

о

<9

Таблица 1

ОТ их ^"■^^(даН/мм2) аа (даН/мм3) -1,5 -0,5 0,5 1.5

0,5 0,25 29,6 28,56 0,46

1,5 22,04 23,03 0,02

2,5 0,005 4,92 4,94

3,5 0,005 0,315 0,215

Таблица 2

50 60 100 300

0,1165-10-ю 0,2.88-10-11 0,217-10-12 0

0,367-10-9 0,710- 10—ю 0,105-10-ю 0,108-10-12

0,19-10-8 0,293.10-е 0,393-10-ю 0,52-Ю“12

0,139-10-7 0,564-10-« 0,62-10-9 0,56-10~11

0,904-10~7 0,316-10—7 0,99-10-8 0,154*10-0

0,278-10-е 0,163-10-6 0,344*10-7 0,500-10-9

силу чего для данной таблицы может быть предложена следующая приближенная формула:

здесь А =5:10 10; вт дано в метрах в секунду, ¿ — в метрах.

Интересно отметить, что число полных циклов для всех исследованных режимов сохраняется приблизительно постоянным и равным -—-100 циклам за 12,56 с.

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

ЛИТЕРАТУРА

1. Слобин Б. 3., Т р о ф и м о в О. Ф. Статистический анализ измерений случайной нагруженности для оценки накопления усталостного повреждения. „Вестник машиностроения“, 1966, № 10.

2. Вентцель Е. С. Теория вероятностей. М., .Наука“, 1969.

3. Трощенко В. Т. Усталость и неупругость металлов. Киев, „Наукова думка“, 1971.

4. Болотин В. В. Долговечность конструкций при квазиста-ционарных режимах напряжений. „Инженерный сборник“, Изд. АН СССР. 1960, т. 29.

5. Б о л о т и н В. В. Об оценке ресурса конструкций при дейст-: вии случайных нагрузок (сб. статей „Расчеты на прочность“, вып. 9).; М., Машгиз, 1963.

6. Райхер В. Л. Гипотеза спектрального суммирования и ее применение для усталостной долговечности при действии случайных нагрузок. Труды ЦАГИ, вып. 1134, 1969.

7. Фын Я. Введение в теорию аэроупругости. М., Гостехтеориз-дат, 1959.

8. Тейлор Д. Нагрузки, действующие на самолет. М., „Машиностроение“, 1971.

9. 3 е л ь д о в и ч Я. Б., М ы к и с А. Д. Элементы прикладной математики. М., „Наука“, 1967.

10. Малинин Н. Н. Прикладная теория пластичности и ползучести. М., „Машиностроение“, 1959.

11. Лейбензон Л. С. Курс теории упругости. М., ОГИЗ, 1947.

12. П и с а р е н к о Г. С., Лебедев А. А. Сопротивление материалов деформированию и разрушению при сложном напряженном состоянии. Киев. „Наукова думка“, 1969.

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