Добавил:
Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз: Предмет: Файл:

книги / Применение метода конечных элементов к расчету конструкций

..pdf
Скачиваний:
8
Добавлен:
12.11.2023
Размер:
20.57 Mб
Скачать

конструкции.

Диссипация энергии в конструкции при ее колебаниях складывается из суммы энергий, поглощенных по каждой фор­ ме собственных колебаний. Форма собственных колебаний, определяющая наибольшее поглощение энергии колебаний, за­ висит от частоты возбуждающей колебания конструкции силы. Демпфирование по одной собственной частоте можно оценить, измеряя в эксперименте декремент логарифмического затуха­ ния для соответствующей формы собственных колебаний. В общем случае необходимо в одном КЭ иметь п значений коэф­ фициентов демпфирования (1~1,2,3,.. .п), если п - число степеней свободы структуры.

Описанный выше способ учета демпфирования является приближенным. Существуют и другие способы, один из кото­ рых заключается в представлении матрицы демпфирования в виде линейной комбинации матриц касс и жесткости [7]:

D=ceM+j3K.

(6.2.7)

Эта модель демпфирования соответствует модели, предложен­ ной Релеем. Константы ос и |3 вычисляются по двум значениям коэффициентов демпфирования, определяемым из эксперимента для двух различных частот собственных колебаний. Считает­ ся, что демпфирование является пропорциональным, т.е.

 

 

 

^kD^l=^uk^k^ kl'

 

(6.2.8)

где

-

собственный

вектор,

соответствующий k-й форме

собственных

колебаний;

-

круговая

частота собственных

колебаний

по k-й форме;

-

коэффициент демпфирования;

$ki

- символ

Кронекера:

 

 

 

 

 

 

 

6

г 0,

если

i*j;

 

 

 

 

=•

если

i=j.

 

 

 

 

kl

[ 1,

 

Подставляя

формулу (6.2.7) в

(6.2.8),

получим:

 

 

 

р£(аМ+/ЗК) ^ = 2 ( 0 ^ 6 ^ .

(6.2.9)

Собственные векторы обладают следующими свойствами:

vkM^k=1' ^ МП = 0 , ^kK^k=wk' HcKvl=0*

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

С учетом этих свойств из равенства (6.2.9) следуют два уравнения с неизвестными а и 0:

а-/3и£=2ыкт7к^

2 с£НЗа)^=2а>^71^.

Из этих уравнений определяем ос и 0:

2

2

/3=

V l ' V k

2 2

wr

“,k

 

Wl“"k

ПРЕОБРАЗОВАНИЕ МАТРИЦ МАСС И ДЕМПФИРОВАНИЯ ПРИ ПЕРЕХО­ ДЕ В ГЛОБАЛЬНУЮ СИСТЕМУ КООРДИНАТ. Матрица масс Mi и мат­

рица демпфирования

вычисляются в местной (локальной)

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

где

- матрица направляющих косинусов углов, вычисляе­

мая по формуле (2.2.5).

б.З. СОБСТВЕННЫЕ КОЛЕБАНИЯ СООРУЖЕНИЙ

Условие равновесия i-ro КЗ имеет вид (формула(2.1.8)):

 

 

Kivi=fi'

 

(6.3.1)

где

- матрица жесткости i-ro КЗ;

- вектор

перемеще­

ния

узлов;

- вектор внешних сил,

определяемый

формулой

(2.1.10).

В случае собственных колебаний в формуле (2.1.10) удельная объемная сила р равна:

p=-pGi

2

(6.3.2)

dt

r

где р -плотность материала; G. - матрица функций формы

i-го КЭ; d

2

v/dt

2

^

 

 

- ускорение узлов.

Все остальные

слагаемые равны нулю. Тогда:

(6.3.3)

или

d V

9

(6.3.4)

dt

где

- матрица масс i-ro КЭ.

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

и V.

(6.3.5)

KV+M ■ 0- =0.

дХг

Уравнение (6.3.5) описывает собственные колебания структуры. Изучение собственных колебаний является важным для исследования поведения структуры придействии на нее внешних нагрузок, зависящих от времени.

Отметим отличие собственных колебаний структуры от свободных. Матричное уравнение свободных колебаний струк­ туры включает слагаемое, соответствующее демпфированию:

dv

d^v

(6.3.6)

Kv+D---+М ---~ =0,

dt

dtz

 

где D - матрица демпфирования; dv/dt - вектор скорости узлов структуры.

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

лизация колебаний структуры в виде уравнения (6.3.5).

• I

Полагая v=p е1а) и подставляя.его в уравнение (6.3.5)

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

К<р-\Н<р,

(6.3.7)

где А=ш Г о) - круговая частота колебаний.

Собственные значения уравнения (6.3.7) расположим в порядке возрастания; OsA^sA^ ...=;АП , где n-порядок сис­

темы (6.3.7).

Каждому собственному значению А^ отвечает собственный

вектор фк. В задачах, связанных с техническими приложе­

ниями, "как правило, вычисляются несколько пар (А^,?ь)

(i— 1/2,...,р).

Причем р<п. Для р пар уравнение (6.3.7) можно предста­ вить в виде:

КФ=МФЛ,

(6.3.8)

где Ф - матрица, столбцы которой есть собственные векто­

ры; А - диагональная матрица собственных значений.

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

Приведем наиболее распространение методы решения урав­ нения (6.3.7).

МЕТОД ВЫЧИСЛЕНИЯ КОРНЕЙ ХАРАКТЕРИСТИЧЕСКОГО ОПРЕДЕЛИ­ ТЕЛЯ. Характеристический определитель для уравнений

(6.3.7) имеет вид:

 

 

Р (А )=det(К-АМ)=0.

(6.3.9)

Алгоритм расчета состоит из 3 этапов.

двух

Этап 1 - * вычисление значений

полинома Р (А ) для

выбранных значений А ^ _ ^ и А ^

(к=1,2,3,...).

 

При выбранном A^j матрица К - А ^ М разлагается на

ниж­

нюю треугольную матрицу и с единичной диагональю и верх­ нюю треугольную матрицу S:

K-A(k)M=u s ,

(6.3.10)

откуда

 

п

 

 

 

(6.3.11)

det(К-А. М )=G S ...

 

1к )

11

 

Этап 2 -вычисление первого собственного значения с

применением метода хорд

 

 

Пусть А (к-1) и А (к)

два

приближения к А1

такие, что

^(k-l)<A(k)^Ai (рис. 6.3.1).

Рис. 6.3.1. Метод хорд для вычисления собственных значений.

По формуле улучшенного метода хорд вычисляем:

х ( к + 1 Г х < к Г г Т7Г

Р ( Х

(

к ) }

Г (Л( к ) ' Л(к-1)>'

(6.3.12)

 

 

 

 

 

 

Е(А( к ) 1 г ,А ( к - 1 ) )

 

 

 

 

 

где гг=2. Пусть

г=2

и X^ +^j<Xa, где

Аа - абсцисса локаль-

ного

минимума

Р(А).

 

отличаготся

в

третьем знаке,

тогда

Если X(к-1) И X( к )

значение г

удваивается.

 

 

 

 

 

 

 

 

Смена знака Р(А) контролируется с использованием свой­

ства

последовательности

Штурма,

которое

заключается в

том,

что

по

изменению

 

знака

Р(А)

между

X(к-1)

и

X(к)

можно разделить собственные значения.

различаются

Переход

к

этапу

3,

 

если

A^+i)

и

только в 7-м знаке.

 

собственного

вектора

методом

об-

Этап 3 - вычисление

ратных итераций.

Обозначим через х^ начальное приближение к собствен­ ному вектору Pj, причем; у1=Мх1>

Для к=1,2,3,.... вычисляем:

Kxk+i=yk'

Ук+1_Мхк+1

_

Хк + Л

7j{xk+i)»

 

/~ т хк+1ук+1

 

yk+l=

ук+1

 

^ * к + Л + 1

 

 

 

При к-*ю

 

 

Обозначим A ^ + I ) . . ^ ^ ^ )

Проверяем при каждом значении к неравенство:

 

|л<к+1>-х<к >|

 

 

 

ТН+ТУ

^ к Г 28,

 

 

 

 

где s-количество верных цифр в числе.

При

выполнении этого

 

неравенства:

 

Xls,,(xk+1>' *1=

к+1

 

^ хк+1ук+1

 

 

 

 

Для

вычисления пары

(\2,<р2) и последующих пар перехо­

дим к

новому полиному:

 

""

 

 

Р ( Л)

 

 

 

 

 

 

 

 

+ ( ^ i >

 

 

 

i=i

1

где

- ранее вычисленные собственные значения (1*1*2»

 

При вычислении

каждой

пары повторяются этапы 1-3.

МЕТОД ИТЕРАЦИЙ В ПОДПРОСТРАНСТВЕ.Рассмотрим вычисление р наименьших собственных значений и векторов уравнения (6.3.8). Метод расчета заключается в следующем:

-вычисление компонент q начальных векторов (q>p)‘;

-выделение "наилучшего" собственного значения и "наилучшего" собственного вектора среди q итерационных векто­ ров методами обратных итераций и Ритца;

-проверка правильности отыскания собственных пар с ис­ пользованием свойств последовательности Штурма после завершения итерационного процесса.

Итерация в подпространстве заключается в следующем.

Пусть X,, (к=1,2,... ,р) - р начальных векторов, обра-

Л .

зующих подпространство Е^. Выполняются обратные итерации:

КХк+1=МХ^,

(6.3.13)

где к*™ 1,2,3,...

р-мерное подпростран-

Тогда р векторов Хк+1 порождают

ство Ек+1.

Ек+1 СХ°ДИТСЯ К Ею/

если начальные векторы неортогональны к Еда.

Вычисляем:

Xk+l_Xk+lSk+l'

(6.3.14)

где Sk+1 - верхняя треугольная матрица, для которой:

(6.3.15)

где I - единичная матрица.

В пределе Xfc+1-*f, Sk+1-»A, которые удовлетворяют уело-

виям:

 

ФТКФ=Л,

(6.3.16)

ФТМФ=1.

(6.3.17)

Последовательность расчета состоит из следующих этапов. Этап 1. Выбираем в качестве первого столбца матрицы

| диагональ матрицы М, в качестве компонент "остальных

столбцов - единицы, соответствующие степеням свободы с наименьшими отношениями kjj/nnjj/

Выбираем q=min{2p,p+8>.

Этап 2. Для к=1,2,3,... подпространство

преобра­

зуется в подпространство Ек+1 по формуле:

кхк+1-мхк.

Определяются проекции операторов К и М на Е^+1:

кк+1=*к+Л + 1 '

“k+i^k+i^k+i'

причем на каждом к+1 шаге решается вспомогательная задача на собственные значения:

Kk+lQk+l=Mk+l®k+l-k+l‘

Улучшаем хк+1: хк+1=хк + А + 1 .

В конце каждого шага цикла производится проверка на окончание цикла.

Этап 3. Проверка правильности нахождения р собственных пар с использованием свойства последовательности Штурма.

 

МЕТОД ЛЛНЦОША• Идея этого метода заключается

в преоб­

разовании исходной матрицы к трехдиагональной

матрице,

для

которой затем решается задача на собственные значе­

ния. Рассмотрим

использование метода Ланцоша для вычисле­

ния

собственных

значений и собственных векторов

матрично­

го уравнения (6.3.7). Умножим обратную матрицу М * на обе части этого уравнения слева .Тогда придем к стандартной проблеме собственных значений:

Aip—Xft

(6.3.18)

где А=М~1К. Так как матрицы М и К - симметричные, то мат­ рица А будет также симметричной.

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

*k+A =A? - V k - i k - A - i '

<6 -З Л 9 >

где к-1,2,3,...п; Ф0=0.

Соотношение (6.3.19) следует из матричного равенства:

ФТ*АФ,

ГД®

•**'-n^ “ °Рт°гональная матрица.

Трехдиагональная матрица Т имеет вид:

 

 

&1

0

0

0

 

«2

^2

0

0

 

Р2

*3

 

0

 

симметрично

 

 

Алгоритм

расчета элементов матрицы т следующий.

Задается

начальный вектор

причем 11^11=1.

Для к=1,2,3,...гп выполняются вычисления:

1) чк-А*к ;

2 > “k= S k ^ k ;

3 > гк=чк-4к“к - * к - А - 1 !

4)Эк- « г к »;

5)-к+1=гк/рк

Если |3^=0 на каком-то шаге, то в качестве Ф^+1 выби­

рается любой единичный вектор, ортогональный к уже вычис­ ленным векторам fj (j=l,2,3,... ,k).

Собственные значения и собственные векторы матрицы Т,

являющиеся собственными значениями и собственными векто­ рами матрицы А, вычисляются следующим образом:

- для каждого к=1,2,3,... решается уравнение:

(6.3.20)

V k * ® k V

- вычисляются аппроксимирующие векторы:

V * kxk'

где Фк=[Ф1гФ2г***/5к ]?

- при

выполнении неравенства

llA^k“0k i V

про

|=s— ------ -—

 

 

 

 

"^к1

 

цесс

вычислений

заканчивается

i_i? i=l,2,3,...).

Вычисленные

и

являются

аппроксимирующими

парами

Релея-Ритца для

собственной пары

(р. ) исходном

матри-

цы А.

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

6.4. ВЫНУЖДЕННЫЕ КОЛЕБАНИЯ СООРУЖЕНИЙ

Записывая в правой части уравнения (6.3.6) вектор си­ лы, зависящей от времени , получим уравнение вынужденных колебаний структуры, моделирующей сооружение:

,dv

 

d2v

(6.4.1)

Kv+D---+М---=f(t),

dt

 

dt:

 

где К, D, M - матрицы

жесткости, демпфирования, масс

структуры; v, dv/dt,

d

v/dt

- векторы перемещения,

скорости и ускорения узлов (производные по времени в

дальнейшем будем обозначать у, v); f(t) - вектор внешней нагрузки.

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

Методы решения уравнения (6.4.1) делятся на две груп­ пы: прямого интегрирования и разложения решения в ряд по собственным формам колебаний структуры.

Считается, что в начальный момент времени t=0 известны векторы перемещений vQ, скоростей vQ, ускорений v Q. Ре­

шение строится для интервала времени от 0 до Т .

В методах прямого интегрирования временной отрезок (0,Т) разбивается на п равных интервалов At: