Динамики упругих конструкций с локальными нелинейностями Текст научной статьи по специальности «Физика»

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

№ 2

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

Г. В. Вронский, В. Д. Ильичев

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

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

— синтез и редуцирование матричных динамических характеристик летательного аппарата без нелинейных элементов (см. [1]);

— формирование отдельных частотно-временных уравнений нелинейных элементов методом исключения;

— формирование дискретного канонического спектрального представления стационарной случайной внешней нагрузки;

— расчет колебаний нелинейных элементов и линейной части конструкции ЛА;

— расчет матричных частотных характеристик конструкции от пространственно-дискретизированной случайной нагрузки с выходом по напряжениям;

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

Методы решения перечисленных задач включают, прежде всего:

— метод синтеза и редуцирования матричных динамических характеристик ЛА [1];

— метод конечных элементов (МКЭ) [2];

— метод быстрого Фурье-преобразования (БПФ) [3, 4];

— метод статистических испытаний (Монте-Карло) [3, 5].

Известно, что МКЭ в задачах динамики приводит к системам

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

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

Рассмотрим Фурье-преобразование уравнений движения дискретной системы в частотной области в виде

бц (/со) С?! (/со) + ЕС12 02 (/со) — Рх (/со); Ю21 & (/ш) + Р{022 [І, Ог (/)]} = Рг (/со);

(1)

здесь (?1 (/<»)•—комплексные амплитуды колебаний всех степеней

свободы дискретной схемы ЛА, исключая степени свободы стыка с нелинейными элементами и степени свободы, отнесенные собственно к нелинейным элементам, составляющие вместе подвек-тор <32;

Г 1^22 Ц, Оа (/)]} — Фурье-преобразование от группы нелинейных уравнений, выражающих движение нелинейных элементов при —

622 К <3:>(^)] — нелинейный дифференциальный оператор;

РД/со), Р2(/со)-—комплексные амплитуды внешних случайных ста' ционарных нагрузок, приложенных соответственно в узлах со степенями свободы О, и Без ограничения общности будем считать связи между большой первой и небольшой второй группами уравнений (1) статическими линейными связями, так как степени свободы стыка в расчетной схеме МКЭ всегда можно связать с основной конструкцией безынерционными упругими конечными элементами (см.

* Или, что еще более эффективно, с использованием микропрограммного специализированного процессора БПФ.

6— Ученые записки № 2 81

фиг. 1, а, б). Тогда 6,2 = 621 — вещественная матрица, что обозначено отсутствием у нее аргумента /со; I — параметр величины связей между и (32. Определив (^(ш) из первой группы уравнений (1)

& (/«>) = С"1 (/со) [- 60,2 0"2 (Ы) 4- Р[ (га.)], (2)

исключим его из второй группы

Г {622 [/, 0_, ^)]} — ?2 Д622 (1т) 0, (/со) = Р2 (ш) — 6Д Р2 (/<»),

где

ДС?22 (/<°) — 621 Си* (/ш) <^?12,

ДМ, (/со) = 621 бц1 (*'“>) £*1 (*’“)•

(3)

(4)

(5)

Л

Фиг.

Применяя к (3) обратное преобразование Фурье -Г-1, получим уравнения движения нелинейных элементов в системе ЛА во временной области

012 [*, $2(*)1 —е2^-1 {АОаз (л*и) /=■ [О;(0]} = /^~1 {Р2(/*«) — Д,(го,)}. (6)

Движение конструкции ЛА определяется теперь в соответствии с (2) соотношением

^ (т) — (га>) — $д (/со), (7)

где С?? (/со) = Оп1 (/со) Р1 (/со) — комплексные амплитуды <3, при

фиксированных степенях свободы нелинейных элементов и стыков _ _ _ (§2 = 0);

Д(3, (/ш) = биЧ*'10) 012 /;l[(Q2 (/)]— приращения амплитуд СМпри ?=1)

за счет движения <32(/).

Зная перемещения С2-> с помощью матричных характеристик 6, вычисляемых в МКЭ, можно рассчитывать напряжения з в расчетных точках различных конечных элементов, так как б связывают перемещения <3 с напряжениями о [1]:

а (/ш) = 6 <3 (/со) ИЛИ а(£) = 0(3(/').

(8)

Статистически разыграем какие-либо из реализаций случайных стационарных внешних воздействий V внешних нагрузок Рх{ш), Р2(г’ш), и, возможно, кинематических случайных воздействий, наличие которых в операторе С22 [/, <32(0] учитывается аргументом Ь. (Пример кинематического воздействия, вызванного неровностями аэродрома, показан на фиг. 1, а).

Представим случайную внешнюю нагрузку, считая ее нормальной, в виде канонического разложения:

Р (№>)

здесь с — случайный вектор с V комплексными независимыми случайными нормальными компонентами, имеющими среднее 0 и дисперсии 1; а(гш) — матричная (неслучайная) каноническая характеристика внешних нагрузок.

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

022 V, & (*)] - Р-1 {021 т (ш) (*)]} = Г-1 {а, (1«) - Ю21 (*»)} 3 =

= Т7-» {а*2 (гш)} с (9)

и

0'1(г<в) = (о (ш)с — 6?(го»)Г [<32(0]> (Ю)

где т(№), ?(*“) являются решениями следующих алгебраических матричных уравнений

Ои(№)т(го>)= 012; Ои (гсо) ср (т) = а1 (ш). (11)

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

2. Общая структура алгоритма частотно-временного метода установления. После формирования и разрешения уравнений (И) „большой" линейной части системы (2) относительно ч(ш) и ?(*’“) задача сводится к интегрированию по времени нелинейной системы уравнений (9) небольшого порядка для заданной реализации внешних воздействий с какими-то начальными условиями.

Если Сп (гоо) не зависит от со (что имеет место, например, если фиксировать узел 1 в системе А, фиг. 1), то (9) принимает, как легко убедиться, вид системы обыкновенных дифференциальных уравнений

о22 к я2 т - оа1 ой1 о12 (о = р* (а (12)

где Р* (0 = Р2 (г) - Ш21 Оп1 Р, (*).

Если 6=0, то „большая" и „малая" подсистемы уравнений не связаны, <$! и С?2 вычисляются независимо друг от друга.

Будем решать (9) численно, комбинируя численное интегрирование по времени и численное дискретное быстрое преобразование Фурье (БПФ) для переходов из временной области в частотную и обратно на каждом шаге (или на каждом /-м шаге) интегрирования по времени. Если БПФ осуществляется на базе N равноотстоящих по времени значений интегрируемой функции, то на первых Мр2 шагах интегрирования положим в (9) 6 = 0. На последующем шаге полагаем 6 = 1, и дискретным преобразованием Фурье N отсчетов функций Ог (£*); к — 1 -+- N получаем спектр С}2(1и>). При этом С}2(£) задается путем экстраполяции, в частности, при £>N/2 можно полагать ^2(^) = 0. Далее, согласно (6), после применения обратного БПФ, получим N отсчетов функции

Д ((К) — 77-1 {(ш) С}., (»}, £ = 1 ч~ N.

где, вообще говоря, Д(^)^0 при &>М2 и, следовательно, может быть сделан следующий шаг интегрирования уравнений (9) методом Рунге — Кутта, при 6 = Iй. Следующие N отсчетов функций (/*) Для применения БПФ берутся в точках & = 2 (/V + 1) и т. д.

Экономичность такого алгоритма может быть повышена за счет экстраполяции Д (/) не на один, а на несколько шагов и, сле-

* Для ускорения процесса установления решения целесообразно переходить от 5 = 0 к ; = 1 не за один, а за несколько шагов интегрирования по времени.

довательно, уменьшения числа БПФ, или за счет упрощения процедуры БПФ для экстраполяции Д(/) на один шаг.

Интегрирование уравнения (9) проводится до „установления" решения, когда возмущения, вызванные начальными условиями и введением в (9) членов с множителем 6, станут пренебрежимо малыми. Установившееся решение (9) позволяет определить искомое движение всей системы с помощью (10). Время установления можно оценить путем сравнения, например, решения (1) при 6 = 1, полученного обычным методом Рунге — Кутта, и решения этого же уравнения в форме (9), полученного частотно-временным методом.

3. Пример расчета. Предлагаемый метод был апробирован на модельной задаче. В качестве конструкции рассматривалась однородная балка А на нелинейном „шасси" В, катящемся с постоянной скоростью V по неровной поверхности (фиг. 1, а).

Движение этой системы описывается дифференциальными уравнениями с частными производными:

EIu™ (t, х) + mu (t, х) + 80 (х — a) {k [и (t, а) — q (£)} + ри (t, а) —

— -ци'х (t, а)} = 0; (13)

Mmq(t)^r k[q(t) — u(t, а)] + Ф(/, q, q) = 0; (14)

здесь u(t,x) — вертикальное перемещение балки;

q (t) — вертикальное перемещение массы Мш (точки 2);

El = EI0 ^1 + i -~j — комплексная жесткость балки;

т — погонная масса балки;

80 — логарифмический декремент колебаний; а — расстояние от конца балки до точки 1\ k — жесткость пружины, соединяющей точки 1 и 2; Р и т) — коэффициенты нестационарных аэродинамических сил от горизонтально расположенной жесткой несущей поверхности, прикрепленной в точке /;

Мт — масса „шасси", которая сосредоточена в точке 2; Ф (t, q, q) — нелинейный оператор, численное значение которого равно силе, возникающей в стойке шасси;

80 (х — а) — дельта-функция.

В рассматриваемом здесь случае полагалось, что

ф(U Я. ?) = *[?(*) — 8(0] + — *(*)] +

+ *н [Я (*) - 3 (013 sign [q (t) - 3 (/)] +

+ са [q (i) — § (012 sign \q (t) — 8 (/)], (15)

где o(^) — функция неровности поверхности (заданная либо полученная методом статистических испытаний из известной спектральной плотности неровностей);

*. *н, ci сп — коэффициенты, учитывающие вклад в значение силы Ф от членов линейной и нелинейной жесткости и линейного и нелинейного демпфирования соответственно.

Представим вертикальные перемещения балки и^,х) в виде разложения по п формам собственных колебаний балки [6]:

и (і, х) = X/; (■*) <7; (О-

(16)

/=і

Подставляя (16) в уравнение (13) и ортогонализируя невязку ко всем формам, получим систему обыкновенных дифференциальных уравнений вида:

Мп §,(0 + (Яп+Яп) 4(*)-МС„ +В?1)^У) + О12дМ = 0-, (17)

здесь Мп — матрица обобщенных масс; (/),, + ^п) — матрица обобщенного конструкционного и аэродинамического демпфирования; (бц 4 Вап) — матрица обобщенных конструкционной и аэродинамической жесткости; С^^) — вектор перемещений [см. (1)]:

Я\ (О

(0 =

Яг (О

Яп (О

(18)

Матрица Оп в уравнении (1) выражается через значения форм в точке х — а

А(а)

/г (а)

0,2 =

к

їп(а)

21 ♦

(19)

Уравнение (14) переписываем в виде [см. (19)]:

Огі Ql (Ї) {Мш д (Ї) кд (і) + Ф (і, д, д)\ — 0. Преобразованием Фурье эта система приводится к виду (1),

(20)

где

Оц (г'ш) =

(И =

і)п) Ч- + Вц,

, <32 (ш) = д (гш),

Ми Ш2 Ш (£)ц

Ч\ (*“)

Яг (П Яп («“)

022 [*, С?2 (0] = (*) 4- кд (0 + ф (/, <?).

Для анализа предлагаемого метода численно находилось решение уравнения рассматриваемого примера двумя способами: методом Рунге — Кутта интегрировались уравнения вида (1), переведенные во временную область; частотно-временным методом установления интегрировались уравнения вида (9).

В расчете учитывалось три тона изгиба балки и ее поступательное вертикальное перемещение. Собственные частоты изгиба балки соответственно равнялись 12,6; ‘34,5; 67,6 рад/с. Задавался гармонический профиль неровностей 8 = агзт2х, вызывающий гармонические колебания. При постоянной скорости пробега V он представляется в виде

0 при Ь < О агэт о>г Ь при / ]> 0; здесь а>8 = 1/2.

При значении V = 100 км/ч исследовалось гармоническое возбуждение при различных амплитудах кинематического воздействия Оо и <»б = 15 рад/с, а также полигармоническое возбуждение с частотами 12, 15 и 18 рад/с.

Результаты расчета перемещений д (/) массы Мш по времени при гармоническом воздействии представлены на фиг. 2. По оси абсцисс отложен номер шага интегрирования, интервал БПФ был равен Т = УУД£, где N=256, М = 0,0082 с.

Сравнение полученных результатов позволяет сделать вывод, что по прошествии некоторого переходного периода, когда наступает режим установившихся колебаний, решения, полученные обоими методами, сближаются. При заданных условиях можно считать, что решения отличаются между собой не более чем на 3% при /2?2Тбпф, где Тбпф — интервал преобразования функции д (£) с помощью быстрого преобразования Фурье (БПФ). При этом порядок линейной части системы практически не влияет на время решения частотно-временным методом.

Для оценки точности и времени установления частотно-временного метода будем сравнивать значения последовательно пронумерованных экстремумов функии <7(^), полученных при гармоническом возбуждении для различных значений параметров нелинейностей шасси частотно-временным методом и методом Рунге— Кутта. Экстремумы функции д(/) можно характеризовать двумя числами: временем Ьк, соответствующим этому экстремуму, и его значением дУк), где 6 — номер экстремума.

Как показали расчеты, значения полученные этими двумя методами, близки между собой и практически совпадают при ^>7’бпф (отличие не более чем на Д/). Поэтому для характеристики точности частотно-временного метода бралось отношение значений <7(?*,), полученных методом Рунге—Кутта и частотно-вре-

' Чр.К и

менным методом--------гг\ в зависимости от номера экстремума к

Ячъ

(фиг. 3).

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

Фиг. З

Фиг. 4

не более чем на 3%, что согласуется с рассмотренным ранее случаем. Число точек N при этом не имеет значения, если частоты, соответствующие пикам амплитудно-частотной характеристики системы и спектра возбуждения, не превышают частоту Найквиста

На фиг. 4 приведена корреляционная диаграмма <7р. к(0 и дчв^) при полигармоническом воздействии на нелинейное шасси (нелинейны жесткость и демпфирование). Интервал применения БПФ соответствовал N=256 отсчетам, цифры на кривой означают номер шага интегрирования. Эффект установления решения выражается в том, что корреляционная диаграмма с возрастанием номера шага стремится к биссектрисе координатного угла.

Анализ полученных результатов показывает, что значения амплитуд колебаний, полученных по методу Рунге —Кутта, превышают по модулю соответствующие значения, полученные частотно-временным методом, приблизительно на 2%. Это различие может быть вызвано рядом причин, например специфическими свойствами дискретного преобразования Фурье на конечном интервале по сравнению с интегралом Фурье, погрешностями в определении частотной характеристики линейной части системы, различием порядков решаемых систем. Для выяснения причины его появления требуется дальнейший анализ.

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

ЛИТЕРАТУРА

1. Ильичев В. Д. Матричные методы синтеза динамических и упругих характеристик линейных неконсервативных конструкций. .Ученые записки ЦАГИ“, т. 6, № 2, 1975.

2. Зенкевич О. Метод конечных элементов в технике. М., „Мир", 1975.

3. Вронский Г. В. Применение быстрого преобразования Фурье и метода Монте-Карло для расчета повторяемости нагрузок и усталостного повреждения элементов авиационных конструкций при колебаниях от стационарной случайной внешней нагрузки. „Ученые записки ЦАГИ“, т. 7, № 4, 1976.

4. Голд Б., Рейдер Ч. Цифровая обработка сигналов. М., „Советское радио", 1973.

5. Бусленко Н. П., Голенко Д. И., Соболев И. М., Срагович В. Г., Шрейдер Ю. А. Метод статических испытаний (метод Монте-Карло). М., Физматгиз, 1962.

6. Б и с п л и н г х о ф ф Р. Л., Эшли X., Халфмен Р. Л. Аэроупругость. М., Изд. иностр. лит., 1958.

Рукопись поступила 22/\г1 1976 г.