Численное исследование динамического поведения вращающихся конструкций Текст научной статьи по специальности «Физика»

УДК 539.3: 534.1

Н.А. ШЕВЕЛЕВ, И.В. ДОМБРОВСКИЙ Пермский государственный технический университет

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ДИНАМИЧЕСКОГО ПОВЕДЕНИЯ ВРАЩАЮЩИХСЯ КОНСТРУКЦИЙ*

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

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

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

Построены амплитудно-частотные характеристики неконсервативной системы для различных значений параметра угловой скорости. Значения последнего соответствуют устойчивому и неустойчивому режиму.

Постановка задачи о свободных колебаниях

В трехмерном евклидовом пространстве упругое тело занимает объем V, ограниченный поверхностью S. На части поверхности Su заданы граничные условия

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

Математическая постановка задачи включает в себя уравнения движения

[1], [4]:

(1 + m)grad (div (u )) + mAu = p[ü + 2щХ u + щх (щх u)], (1)

где u - вектор перемещений во вращающейся цилиндрической системе координат r, z, фє V, ось z которой направлена по вектору щ; p - плотность материала тела;

1, m - параметры Ламе.

* Работа выполнена при финансовой поддержке РФФИ проект № 09-08-99121-р_урал_офи.

© Шевелев Н.А., Домбровский И.В., 2009

Граничные условия в перемещениях на части поверхности Еи и в напряжениях на части поверхности Еа ( Е = Е и + Еа )

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

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

Здесь 0 - первый инвариант тензора деформаций, Е - единичный тензор второго ранга. Для компонент тензора деформаций и вектора смещений имеют место соотношения Коши:

Искомый вектор перемещений для рассматриваемой задачи представляется в виде

где р - комплексное собственное значение, имеющее смысл собственной частоты колебаний или параметра устойчивости.

После подстановки выражения (5) в уравнения (1) и граничные условия (2) получим формулу для отыскания вектора и, который представляет собой форму колебаний, если р - мнимое число, или форму потери устойчивости, если действительная часть комплексного р или действительное р больше нуля, краевую задачу

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

и(г, г, ф, і) = 0, г, г, фє Еи , н а (и ) = 0 , г, г, фє Еа,

(2)

а = 10Е + 2|іє.

(3)

1 Г г є = — (Уи) + Уи

21

(4)

и (г, г, ф, і) = и(г, г, ф)ер,

(5)

(1 + т)Бгаё(ё1у(и)) + дДи = р р2и + 2рщхи + щх(щхи) , г,г,фє V, и(г, 2, ф) = 0, и(г, 2, ф) = 0, г, г, фє Еи,

(6)

8Ла + 84 + 8ЛС + 8Лк = 0,

8Ла = -1 а • 8єс^, 8Л = -рр2 { и • 8и^, 8ЛС =-| р [щх(щх и)] 8и^, 8Лк =-2рр| (щх и )8и^

Здесь 8Ла - работа внутренних напряжений, определяемая возможными

перемещениями, 8Д. - работа сил инерции, 8ЛС - работа центробежных сил, связанных с полем перемещений и, 8Лк - слагаемое, определенное силами Кориолиса. Интегрирование ведется по объему V.

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

Численная реализация

Геометрическая симметрия рассматриваемых тел относительно оси вращения позволяет записать искомый вектор и(г, г, ф) в виде разложения по окружной координате ф:

и =

N

2 п'1 (г, I) С0Б(/ф) + П^ (г, і) 8Іп(/ф)

/=1

N

2 (г, і) С0Б(/ф) + Н’/а (г, 2) БІп(/ф)

/=1

N

2 (г,2) соз(/Ф)+уі(г, 2) ^п(ф)

/=1

(8)

где п, н, V - компоненты вектора и по координатам г, 2, ф соответственно, индекс ^ соответствует симметричным и а - антисимметричным составляющим вектора относительно радиуса ф = 0, суммирование ведется от / = 1 до / = N.

В символической форме вектор перемещений запишется в виде

и = и, + иа . (9)

Если в вариационном уравнении не учитываются слагаемые 8ЛС и 8Лк, связанные с вращением ( ю = 0), то задача определения и распадается на две самостоятельные для и иа, а соответствующие собственным векторам собственные значения р5 и ра

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

На этапе численного построения векторов и и иа применяется процедура метода конечных элементов.

Конечно-элементный аналог вариационного уравнения (7) запишем в виде

(К + рр 2М + рм>2¥с + 2ырр¥к) 8 = 0,

(10)

К = | Вт ШdV, М = | N NdV,

Ж = / N Г^У, Жк = | N {„ЫУ,

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

д 2 д 2

Г0 = = О^у, (11)

0 д^2 дt2

и положить О = ю, то можно уравнение (10) переписать в виде

(К0 +Д,М0 +РоЖ; )8 = 0, (12)

где матрицы с нулевым индексом отличаются от соответствующих матриц из (10) множителями, а р = р0О. Такая замена переменной t и выбор О приводит к нормировке численного значения элементов матриц К, М, Жс, Жк и их равноценному вкладу в результирующие коэффициенты системы (12). Комплексные собственные значения р0 , входящие в систему линейных алгебраических уравнений как

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

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

Разложение включает линейные комбинации конечного числа первых собственных форм колебаний соответствующей упругой задачи без учета вращения [3]:

т

3 = Е Чкчк , (13)

к=1

(чк - собственные формы колебаний). Для их определения решается алгебраическая задача

(К -р12 М) чк = 0, (14)

где 1к - собственная частота колебаний консервативной системы, чк - собственная форма колебаний, соответствующая 1к. При учете (13) система (12) преобразуется в систему

(К + р02М' + р0Гк) д = 0. (15)

Матрицы в (15) имеют размер 2т х 2т, их элементы определяются соотношениями:

К = чК0ч 1, ы9 = чМ0ч■, = ч[Ж,»чу, /,} = 1,2,...т (16)

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

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

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

Примеры расчета

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

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

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

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

нуль. При дальнейшем увеличении ю величина р1, например, становится

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

Рис. Изменение собственных значений в зависимости от параметра угловой скорости

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

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

Вынужденные колебания с учетом вращения и диссипации

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

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

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

U (r, z, j, t) = (U0 + Uc) cos (Wkt) + (U0 + U*) sin (Wkt), (17)

где и0 (г, г, ф) и (г, г, ф) - амплитуды заданного кинематического возмущения;

0.к - частота кинематического возмущения.

Согласно полуаналитическому методу конечных элементов раскладываем искомые перемещения в ряд по угловой координате:

U

U + Uо V + Vo

W + W

Е [(Ucc + U0cc) cos(/j) + (Uac + U0ac) sin(lj)] cos(Wkt)

i=i

¿ [(Vcc + V0c )cos(lj)+(Va + Vac )sin(lj)]cos(Qkt) +

l=1

¿ [(Wcc + W0cc) sin(lj) + (Wac + W0ac) cos(lj)] cos(Wkt)

l=1

Е [(Uc0 + U00) cos(lj) + (Ua0 + Ua0) sin(lj)] sin(Wkt)

l=1

+ ¿ [(Vc0 + VqS ) cos(lj) + (Va0 + V0ao) sin(lj)] sin(Wkt)

l=1

¿ [(Wc0 + W0co) sin(lj) + (Wa0 + W0ao) cos(lj)] sin(Wkt)

l=1

(18)

где и0, Ж0, У0 - компоненты заданного кинематического возмущения.

Будем разыскивать неизвестные симметричные Цсс, Же, Vее, Ц“, Ж, У“ и антисимметричные С/20, Ж°е, У°е, Цх?, Ж™, У* перемещения методом конечных элементов. Используя конечно-элементную модель, запишем для произвольного элемента соотношения (18) в виде

и = Е N (дг )с°8 ) + N (дг )81" 8Ш(^) , (19)

/=1

где (д/ )С°8 =(д/)е+(д0 )е, (д/Г =(д/Г+(д/0)'; (д/ )с, (д/Г - неизвестные узловые

перемещения; (дг0) , (д0) - заданные узловые перемещения.

В отличие от упругого тела в этом случае матица Б имеет вид

D = D* • G* + D0 • K

(20)

С учетом равенства (20) вариационное уравнение в конечно-элементной форме примет для произвольного элемента вид

( t \

К*о

с°8(Йкґ) - | Я(ї -х)с°8(Йкх)ёх (дг )С°8 + К*о (8Іп(Йкґ) -

V -¥ У

| Я(ї - х) 8Іп(Йкх)ёх (дг )8і" + К • К0 с°8(й,ґ) (д1 )С°8 + К • К0 8Іп(Йкґ) (дг )8іп

-Ц2РМ С°8(Ц^) (д/)С°8 -Ц2РМ 8Іп(ЦкҐ) (д/)8ІП -Ю2рРс С°8(ЙкҐ) (д/)С°8 --Ш2рГс 8Іп(ЦкҐ) (д / )”” - 2тЙ,рРк 8Іп(ЙкҐ) (д / )8ІП +

+2тЦРР с°8(Цкі) (д / )С°8 = Ко (с°8(й,ґ)(1 - Г) - Г, 8т(Й^)) (д / Г + (21)

+К*о (8Іп(Цкґ)(1 - Гс) - Г, с°8(й,0)(д/ )8іп +

+К • К0 с°8(Цкґ) (д/ )С°8 + К • К0 8Іп(Цкґ) (д/ )8іп - Ц2рМ с°8(Цкґ) (д/ )С°8 -

-ц2рм 8іп(Цкґ) (д / )8іп - «2рР с°8(Цкґ) (д/ )с°8 - т2рР 8Іп(Цкґ) (д / )8іп --2тЦкррк 8Іп(ЦкҐ) (д/ )8іп + 2тЙкррк С°8(ЦкҐ) (д / )С°8 = 0

где

,-Р(^-х)

\1-а

е

Щ-х) = А-----—, (22)

(ґ -х)1

/• е

Гс = А Г—:---------С°8 Й,0^0,

с J 01-а к ’

Г, = А Г —— 8Іп Й, 0ё0,

* J 01-“ к

К* = Г (В— )Г Ъ*Вееё¥е

Vе

К = Г (В— )Г Ъ0Ве^е.

Т7-Є

Система уравнений (21) распадается на две подсистемы:

(о • К*(1 - Гс)+К • К0 - Й;рМ - т2рр ) (д/ )С°8 +

+(о • К*(1 - Гс)+2трЦж, )(д/ )8І" = 0,

(о • К*(1 - Гс)+К • К0 - П^рМ - т2рр ) (д/ )8Іп -- (о • К'(1 - Гс) + 2трЦЖк) (д/ )С°8 = 0,

(23)

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

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Трояновский И.Е., Шардаков И.Н., Шевелев Н.А. Проблема собственных значений и форм вращающихся деформируемых конструкций // Прикладная математика и механика. - 1991. - Т. 55. - Вып. 5. - С. 857-864.

2. Болотин В.В. Неконсервативные задачи упругой устойчивости. - М.: Физматгиз, 1961. - 339 с.

3. Меркин Д.Р. Введение в теорию устойчивости движения. - М.: Наука, 1971. -

312 с.

4. Колтунов М.А., Кравчук А.С., Майборода В.П. Прикладная механика деформируемого твердого тела. - М.: Высшая школа, 1983. - 351 с.

Получено 1.05.09