Исследование переходных процессов при динамическом нагружении упругих конструкций Текст научной статьи по специальности «Строительство и архитектура»

ДИНАМИКА, ПРОЧНОСТЬ, НАДЕЖНОСТЬ \

УДК 539.3

Ю. И. Клюев, О. А. Нехаевская

ИССЛЕДОВАНИЕ ПЕРЕХОДНЫХ ПРОЦЕССОВ ПРИ ДИНАМИЧЕСКОМ НАГРУЖЕНИИ УПРУГИХ КОНСТРУКЦИЙ

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

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

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

В настоящей статье предложен матричный метод, позволяющий успешно реализовать оба этапа решения. Математические основы метода изложены в работе [1].

Решение линеаризованных уравнений движения. Линеаризованные уравнения движения представляются в канонической форме:

д ( д2 \

- у(м) = (Л + В—^ у(М) + д(м), (1)

где 8 — пространственная координата; t — время; Y(s,t) = = [(и(8,£))т (Р(8,£))т]т — вектор состояния сечения, включающий в себя обобщенные перемещения и соответствующие им внутренние силовые факторы; Л и В — матрицы жесткостных и инерционных характеристик; д — вектор внешних распределенных нагрузок.

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

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

Учитывая, что уравнения теории оболочек являются "жесткими" [1], интервал интегрирования по пространственной координате разбивается на определенное количество участков. Длина участка выбирается из тех же соображений, что и в методе С.К. Годунова.

Удерживая в уравнении (1) только линейные составляющие, содержащие операторы (что справедливо при малой длине участка),

получим для участка вг € [в0, в^] (в дальнейшем индекс г опустим) следующее:

Y(si,t) = = (K (A) + KS0 (A)

K

-1

д2

• B • Ko (A)dedt2 )Y(so,t)+

si

+ KSO(A) / Kn(A) • Q(ß,t)dß. (2)

1

so

Отметим, что решение на участке вг вычисляется независимо от предыдущего участка, поэтому операции ортонормирования (как в методе С.К. Годунова) не требуется.

Для получения корректных решений уравнений с переменными по координате в коэффициентами участок в € [в0, в;] разбивается на М интервалов, внутри интервала коэффициенты считаются постоянными и равными их значениям при вт (вт — координата середины т-го интервала).

В итоге, для уравнений с переменными коэффициентами имеем

m=1

KS0 (A)= П KS:_! (A(em)), m = 1, 2,..., M.

(3)

m=M

Второе слагаемое в выражении (2) вычисляется по следующей формуле [1]:

T = Km • K

M1

• K2 • Bi + Km • K

M1

K3 • B2 • Ki + ... +

o

+ km • km-i • ... • Km+i • Bm • Km-i • ... • К2 • Ki + ... +

+ Км • Вм-1 • Км-2 • ... • К2 • К 1 + Вм • Км-1 • . . . • К 2 • К l, (4)

где

В — Ле В I (Лвт) (В А + А В ) +

Вт ЛетВт + ^ (ВтАт + Ат Вт) +

(Ле )3

+ 3т (ВтАт + АтВт + АтВт) + ... ;

Лет ет ет-1; Кт 1 (А(вт)); Ат А(вт); Вт В(вт)?

т — 1, 2, ...М.

Третье слагаемое в уравнении (2) (вектор частного решения) также вычисляется с помощью нормированного решения однородного уравнения

Ь — Км • Км-1 • ... • К2 • 01 + Км • Км-1 • ... • Кз • 02 +... +

+ KM • KM-1 • ... • Km+1 • Q

+ ... + KM • QM- 1 + Qm , (5)

где

Am (^sm)' Am (Asm)

Qm = (EASm + ^ ^ V~ ^ + m ^ ^ + ...| Q (вт).

Для уравнений с постоянными по пространственной координате 8 коэффициентами а^ матрица нормированных решений уравнения (2) определяется формулой

К0 (А) — Е + ^^ + А'2(8,~ 80)2 +... — £ А (8'7 80 , (6)

• • ¿=0

где Е — единичная матрица; ] — количество удерживаемых членов матричного ряда.

В соотношениях (3)-(6) аргумент £ опущен.

Обозначив матрицу нормированных решений (3) через И, представив И, Т и вектор Ь в блочном виде, решение по координате 8 для ¿-го участка представим в виде

и(8,,£)— (Иц +Т11 Ц2^ и(во,*)+ (и!2+Т12Р(во,*)+М*);

Р(М) — ^21+Т21 и(во,£)^И22+Т22^^ Р(*0, (£).

(7)

Преобразуя соотношения (7) к виду

P(so,t)= (kii+mn) U(so,t)+ (k12+m12¿Л U(s,t)+fi(t),

dt2

'dt2

Р(М)= ( к21+т21 к22+т22И^,*^(¿)

(8)

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

2

M—U + KU = F,

dt2

(9)

где и = [ИЦ и2 ... иП]т, п — число узлов.

Отметим, что матрицы М и К в уравнении (9) получены в результате решения системы дифференциальных уравнений, а не в результате априорного задания функций перемещений (МКЭ).

Они представлены в аналитическом виде — в виде сходящихся матричных рядов. Сходимость таких матричных рядов показана в работе [2].

Запишем уравнение (9) в виде

и + си = Р(£), (10)

где С = М-1К, Р(£) = М-1Е(£).

Интегрирование системы дифференциальных уравнений движения по времени. Уравнение (10) сводится к уравнению первого порядка

¿(¿) = Б • 2(£)С(£), (11)

где Z(t) =

U(t)

U(t)

D

; G(t) =

о

P(t)

0 Е -С 0

Будем считать, что жесткостные и массовые характеристики не зависят от времени.

Решение (11) на интервале времени [0,£] находится так же, как и при интегрировании по координате 8 [1]:

Z(t) = K0(D)Z(0) + K0(D)| [K0(D)]-1 • G(t)dr =

t

= K0(D)Z(0) + у KT(D) • G(t)dT,

(12)

где

Dt D2t2 D3t

\3+3

K0 (D) = E + 1Г + ^т +

3!

+

= e

Dt

(13)

t

KT (D) = E +

D(t — t ) D2(t — t ) D3(t — t )

+

+

+ ... .

(14)

1! 2! 3!

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

Наиболее рациональными с точки зрения времени и точности счета являются вычисления этой матрицы с помощью интеграла Воль-терра [2]

K0(D) = (E + At • D)k, At =

k

Например, при k = 1024, необходимо всего 10 операций умножения матриц.

Рассмотрим использование предложенного алгоритма при ступенчатом приложении нагрузки

G(t) = [ 0 Рт ]т = const.

Используя выражения (12)-(14), получим

U(t)

U(t)

hi h2 Ьб he

Uo U o

+

h3 h4

P,

(15)

где h ... h6 — матричные ряды:

hi = E — 2! Ct2 + 4! C2t4 — 6! C3t6 + 8C418 — .

h2 = Et — 1 Ct3 + j C2t5 — 1 C3t7 + ...; Ьз = Et2 — i Ct4 + 6 C2t6 — 8 C3t8 + ...;

(16)

Ь4 = Ь2;

Ь5 = -с* + 3 С2*3 - 1С3*5 + С4*7 - ...; Ье = Ьь

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

При наличии условий закрепления конструкции для рассмотренного вида нагружения нет необходимости находить частное решение уравнения (11), так как Ь3 = С-1(Е — Ь^, Ь4 = Ь2.

На основе разработанной теории функций от матриц Ф.Р. Гантмахер предложил решение уравнения и + Си = Р(*) в следующем виде:

3

• • ?

-1

U(t) = cos ^v^t] Uo + (л/о) ' sin (Vet) U+

+ (VC

1

sin

VC(t - T)! P(t)dr, (17)

где синус и косинус от матрицы определяются матричными рядами:

н2Р+1 х H2p

sin(H) = > : (-1)^7^^; cos(H) = ^ (-1)p—.

p=0

(18)

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

hi = cos ÍVC^ ; h2 = ÍVc) - sin ÍVC^ ;

ha = (VCj V sin

VC(t - T)

dT.

Аналитические значения интеграла частного решения для нагрузок различного вида. Авторами получены аналитические значения интеграла частного решения (17) для нагрузок широкого класса, в частности:

a/C) / sin

VC(t - T)

PodT = C

1

E - cos (VCt)

Po

— для постоянной нагрузки P(t) = P0;

a/O / sin

VC(t - T)

PiTdT = C

1

Et — ( VC) sin (v^t

Pi

— для линейной нагрузки P(t) = P1t;

а/О / sin

- T) P0 sin(wT)dT =

= (C - Ew2)

21

E sin(wt) - w VC) sin л/Ct)

Po

— для гармонической нагрузки P(t) = Р0 sin(шt).

Решение уравнений движения с учетом демпфирования. Авторами обобщено решение Ф.Р. Гантмахера для исследования колебаний конструкций с учетом демпфирования.

t

t

t

t

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

и + пСИ + СИ = Р(*), (19)

где п =

6

2п2/!

; 0,001 ^ 6 ^ 0,02; 6 — логарифмический коэффициент

затухания; /1 — первая собственная частота конструкции в Гц. Решение (19) при ступенчатой нагрузке Р(*) = Р0 имеет вид

И(*) = ехр(-Ш) [еов(р1*) + Кр-1 б1п(р1*)] Ио + р-1 в1п(р1*)ТИо+ + р-1 (К2 + р2)-1 [р1 - ехр(-Ш)(р1 еоБ(р1*) + N ^(р^)] Ро,

1 __^ ир

где N = -пС; р1 = у/С - К2; ехр(Н) = ^ ~т.

2 р=0 р

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

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

д = 106— (рис.1). Исходные данные: £ = £2 = 6,6708 • 109 Па;

м2

^23 = &13 = ^12 = 2,06 • 109 Па; ^2 = = 0,32; р = 2640 Я = 1м; Н = 0,1м.

На начальном этапе (тт. 1 и 2 (рис. 2)) сфера деформируется с преобладанием мембранных деформаций всюду, за исключением узкой краевой зоны; центральная часть оболочки "еще не знает" о наличии контурного защемления. Однако от контура начинает распространяться изгибная волна, сходящаяся к полюсу. Фокусируясь в центре, она увеличивает деформацию оболочки.

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

кг

мз

Рис. 1. Схема приложения ступенчатой нагрузки к защемленной полусфере

Рис. 2. Осевые перемещения в полюсе полусферической оболочки

Рис. 3. Эпюра продольных перемещений по длине полусферы

чение 2 (в данной задаче ^ 5,5). Результаты расчета совпали с результатами А.В.Кармишина [3]. На рис.3 приведен график продольных перемещений по длине полусферы в момент времени т. 3 на рис. 2.

2. Колебания простейшей модели космического аппарата при приложении ступенчатой нагрузки (рис. 4). Исходные данные: Нх = 0,2 м; Н2 = 0,25 м; Ях = 0,5 м; Я2 = 2 м; М2 = 2000 кг; рама: тонкостенные балки с параметрами: = 0,046 м; ¿2 = 0,042 м; Е = 0,71 х

кг

х 1011 Па; в = 0,271 • 1011 Па; V = 0,32; р = 2780—; Ро = (Мх+М2)8д;

м3

0,02 [Г ] П =2П2Д [ГЦ]-

Рис. 4. Модель космического аппарата при приложении ступенчатой нагрузки

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

Колебания космического летательного аппарата (КА) (рис. 5) при релейном законе управления Рис.6. Релейный закон управления по крену (рис. 6). Первая и вто- по крену рая собственные формы колебаний

КА представлены на рис. 7, а, б. Исходные данные: Мсб. = 56 кг; Мтв.т. = 360 кг.

Рассматривается задача угловой стабилизации космического аппарата с использованием принципа релейных систем (см. рис. 7). В ре-

б

Рис. 7. Собственные формы колебания космического летательного аппарата с солнечными батареями:

первая (а) и вторая (б) упругие формы

Рис. 8. Поперечные перемещения КА с солнечными батареями с учетом (кривая 2) и без учета (кривая 3) упругости солнечных батарей:

1 — управляющее воздействие — момент

лейных системах реактивные сопла функционируют в течение небольших промежутков времени, необходимых для устранения отклонения аппарата от заданного положения. Установившийся режим работы релейных систем — автоколебания. Обычно задача колебаний КА решается без учета упругости солнечных батарей, т.е. аппарат и солнечные батареи представляют собой твердые тела. Авторы исследовали процесс колебаний КА с учетом упругости солнечных батарей (рис. 8).

СПИСОК ЛИТЕРАТУРЫ

1. Виноградов Ю. И., Клюев Ю. И., Образцов И. Ф. Метод решения краевых задач механики деформирования тонкостенных конструкций // МТТ. - 2001. - № 1. - C. 159-166.

2. Гантмахер Ф. Р. Теория матриц. - М.: Наука, 1967. - 575 с.

3. М е т о д ы динамических расчетов и колебаний тонкостенных конструкций / А.В. Кармишин и др. - М.: Машиностроение, 1990. - 288 с.

Статья поступила в редакцию 22.11.2006

Юрий Иванович Клюев родился в 1939 г., окончил МВТУ им. Н.Э. Баумана в 1962 г. Д-р. физ.-мат. наук, профессор кафедры СМ-2 МГТУ им. Н.Э. Баумана. Автор 69 научных работ в области механики деформирования упругих тел.

Yu.I. Klyuev (b. 1939) graduated from the Bauman Moscow Higher Technical School in 1962. D. Sc. (Phys.-Math.), professor of "Spacecrafts and Boosters" department of the Bauman Moscow State Technical University. Author of 69 publications in the field of mechanics of deforming of elastic bodies.

О.А.Нехаевская родилась в 1983 г., окончила МГТУ им. Н.Э. Баумана в 2005 г. Аспирант кафедры СМ-2 МГТУ им. Н.Э. Баумана. Специализируется в области механики деформирования упругих тел.

O.A. Nekhaevskaya (b. 1983) graduated from the Bauman Moscow State Technical University in 2005. Post-graduate of "Spacecrafts and Boosters" department of the Bauman Moscow State Technical University. Specializes in the field of mechanics of deforming of elastic bodies.