Особенности дифференциальной и вариационно-разностной формулировок задачи продольно-поперечного изгиба стержня от сил инерции

Полный текст

Совершенство конструкций техники космического назначения с позиции механики деформируемых сред связано с уменьшением массы, обеспечением прочности и жесткости в связи с активными воздействиями управления и космического пространства. В качестве конструкций, создающих ускорения в наземных условиях, можно назвать центрифугу для тренировки космонавтов, создающую центростремительное ускорение порядка 30 g [1]; многофункциональные высокоскоростные центрифуги для решения медикобиологических проблем [2], развивающие скорость вращения до 500 об/с; в атомной промышленности [3] известные центрифуги имели скорость 2000 об/с. Основоположниками расчета устойчивости конструкций являются Я. Бернулли, Л. Эйлер, Ж. Л. Лагранж, С.Д. Пуассон [4]. По проблемам устойчивости пластин, оболочек, стержневых систем назовем имена С.П. Тимошенко, А.В. Александрова, И.А. Биргера, A. С. Вольмира, Г. С. Писаренко, В.И. Феодосьева, П. М. Варвака, Я.Г. Пановко, Е.П. Попова, Ю.В. Захарова, В.А. Светлицкого, А.Ф. Смирнова, В.В. Кабанова, Л.П. Железнова, В.В. Новожилова, В.З. Власова, Х. М. Муштари, Э.И. Григолюка, В.В. Болотина, B. И. Мяченкова, Ю.В. Немировского, Н.П. Абов-ского, Я. М. Григоренко, Н.В. Пустового, А.Н. Андреева, Л.И. Шкутина, Н.А. Алфутова, С.Н. Кана. Литература по расчету устойчивости стержней обширна. Ее обзор, начиная с книг [5; 6], до работ более современных исследователей, из которых назовем [7; 8], показал, что для защемленных по торцам стержней на действие изменяющихся инерционных осевых нагрузок аналитические решения неизвестны. Поэтому разработанную вариационно-разностную методику численного анализа устойчивости стержней от инерционных воздействий можно считать актуальной. Рассмотрим дифференциальную формулировку задачи. При продольно-поперечном изгибе известно уравнение равновесия в системе координат Oxyz [9] EJwxxxx - 4z = -N ( x) wxx , (1) где E - модуль Юнга материала стержня; J - осевой момент инерции поперечного сечения, ориентированного в осях Oyz; w = w(x) - функция прогиба стержня по направлению минимальной жесткости; N ( x) - функция продольного внутреннего усилия. При решении задачи на собственные значения поперечная распределенная нагрузка qz = 0 . Выполним численный расчет продольнопоперечного изгиба стержня (рис. 1) методом конечных разностей. Для этого нанесем на область стержня длиной l (рис. 2) конечно-разностную сетку 1 < i < n + 3 с шагом Л = l / n . Здесь n - число шагов сетки в области стержня (число целое). Тогда i = 1 и i = i + 3 - законтурные узлы, а i = 2 и i = n + 2 -номера контурных узлов. Дифференциальные операторы в (1) для дискретной области стержня заменим конечно-разностными аналогами [10], что приводит к системе уравнений wi-2 - 4wi-1 + 6wi - 4wi+1 + wi +2 = = -Ni Л2 -2wi + wi+!)/(EJ), (2) 1 < i < n + 3. Если прогибы и углы поворота торцов балки равны нулю w(0) = 0 , w(l) = 0 и 0(0) = dw(x) / dx|x 0 = 0, S(l) = dw(x) / dx|x l = 0 , (3) тогда в (2) принимается 3 < i < n +1 (4) и добавляются граничные условия: w(0) = w2 = 0 , w(l) = wn+2 = 0 , w1 = w3 , wn+3 = wn+1 . (5) Система тождеств (2) представляет обобщенную проблему собственных значений [11] [ A]{w} = s [ B ]{w}, (6) где [A] - матрица жесткости стержня; матрицу [B] назовем матрицей внутренних инерционных сил; {w} = {w3 w4 w5 ,..^ wn-1 wn wn +1} - собст венный вектор; s - набор собственных чисел для 132 Математика, механика, информатика соответствующих форм потери устойчивости. Первой форме потери устойчивости соответствует минимальное собственное число s. Поясним прием (6) на модели балки с редкой сеткой n = 8. В этом случае зависимости (2)-(5) примут вид 0' 0 0 0 1 -4 7 ' EJ " 7 -4 1 0 0 0 4 6 -4 1 0 0 1 -4 6 -4 0 0 1 -4 6 -4 1 0 0 1 -4 5 -4 0 0 0 1 -4 6 0 0 0 0 -4 2 N3 - N3 0 0 0 - N4 2N 4 - N 4 0 0 0 - N5 2N5 - N5 0 0 0 - N6 2 N6 - N6 0 0 0 - N7 2 N7 0 0 0 0 - N8 0 0 0 0 0 0 0 0 0 - N7 2 N8 0 0 0 0 0 - N8 (7) -N9 2 N9 Отметим особенность (7). Здесь матрица B в общем виде не является симметричной матрицей относительно её главной диагонали. Можно переформулировать прием вычисления (6) к виду [B]H = s* [a]{w}, (8) тогда первой форме потери устойчивости будет соответствовать максимальное собственное число s; причем s = 1/ s . Введем параметр qx , имеющий смысл равномерно распределенной нагрузки, приложенной вдоль оси стержня (рис. 1, а). Тогда функция продольной силы (стержень является однажды статически неопределимой системой) равна N ( x) = qx (l /2 - x ) , 0 < x < l. (9) График (эпюру) функции (9) изобразим на рис. 1, б. Подставив (9) в (1), получим дифференциальное уравнение EJwxxxx =-qx (l/2-Х) Wxx , 0 < X < l , которое дает уравнения (7) в следующем виде (10) EJ ' 7 -4 1 0 0 0 0 " W3 ' -4 6 -4 1 0 0 0 w4 1 -4 6 -4 1 0 0 w5 0 1 -4 6 -4 1 0 ■ w6 0 0 1 -4 6 -4 1 w7 0 0 0 1 -4 6 -4 w8 0 0 0 0 1 -4 7 w9 _ qj ' 8X2 " 6 -3 0 0 0 0 0 w3 -2 4 -2 0 0 0 0 w4 0 -1 2 -1 0 0 0 w5 0 0 0 0 0 0 0 ■ w6 0 0 0 1 -2 1 0 w7 0 0 0 0 2 -4 2 w8 0 0 0 0 0 3 -6 w9 (11) Действительно, в (11) матрица внутренних инерционных сил B оказалась несимметричной относительно главной диагонали, также выродилось одно уравнение, что является еще одной особенностью дифференциальной формулировки задачи. Можно уйти от нулевого коэффициента на главной диагонали, добавив или убрав один узел конечно-разностной сетки или подобрав ненулевую нагрузку в узле. Однако несимметричность матрицы внутренних сил остается. Собственные значения и собственные векторы в (11) не могут быть вычислены, даже если задача будет приведена к виду (8). qx qxl qxl 2 б Рис. 1. Стержень, жестко закрепленный по краям: а - нагрузка, действующая на стержень; б - эпюра продольной силы Для классической задачи устойчивости Эйлера с действующей на торец стержня силой P дифференциальная формулировка (1)-(2) проблемы собственных чисел дает решение. Внутренние силы: Ni =-P, i = 3 - 9; матрица B в (6) становится симметричной, причем s = PX2 / EJ = Pl2/(64EJ). (12) Решение (6) даёт минимальное собственное число s = 0,586, подставив которое в (12), получаем минимальное значение силы P (значение критической силы), равное P = 64 EJ s /12 = 37,49EJ /12. (13) Представим возможную форму потери устойчивости в виде функции w( x ) = (1 - cos2nx /1 ) / 2 на интервале 0 < x < l. Подставив производные этой функции l а 4 133 Вестник СибГАУ. 2014. N 3(55) в уравнение (1), приняв N(x) = -P , получаем значение критической силы P = 4EJП /12 = = 39,48EJ /12, практически не отличающееся от результата (13). Для уточнения величины P в численном расчете, следует сгустить конечно-разностную сетку. Заключаем, что при дифференциальной формулировке задачи об ускорении, при переменном внутреннем усилии N(x), решения нет. Переформулируем краевую задачу (1) в интегральную формулировку [12; 13] в перемещениях. Левую часть (1), последовательным интегрированием по частям переведем в функционал Лагранжа: l 2 l Эл (w) = -1 / 2 j EJ(x) (wxx ) dx + j qz w(x) dx +Qz (l)w(l) -0 0 -Qz(0)w(0) -Mx (l) S(l) + Mx (0) 0(0), (14) где 0 < x < l; dx - длина бесконечно малого элемента; Q (l ), Q(0), Mx (l ), Mx (0) - поперечные силы и изгибающие моменты на торцах балки. Преобразуем в (1) правую часть, которая представляет в уравнении равновесия проекцию внутренней силы на вертикальную ось. Запишем выражение возможной работы этой силы на возможном перемещении 5w( x). Дальнейшее интегрирование по частям дает x=l x=l j N(x)wxx5w(x) dx = - j N(x)wx5wxdx + x=0 x=0 +N (l)&(l)8w(l) - N (0) S(0)8w(0). (15) В (15) усилие N ( x ) не варьируем, потому как оно является известной функцией. Присоединив выражение (15) к вариации функционала (14), получаем x=l x=l j EJ(x)wxx§wxxdx - j qz §w(x) dx - Qz (l)Sw(l) + x=0 x=0 +Qz (0)5w(0) + Mx (l ) 5S(l ) - Mx (0) 5S(0) = x=l = j N ( x)wx5wx dx - N (l)S(l)5w(l) - N (0)S(0)5w(0), (16) x=0 добавив к которому главные граничные условия задачи (3), имеем вариационное уравнение l x=l j EJ (X)wxx5wxxdX = j N(X)wx5wxdX . (17) 0 x=0 Сформируем коэффициенты матриц А и В в (6) выполнив варьирование (17), применив операторы варьирования 51 и 52 : l x=l j EJ(x)51wxx52wxxdx = j N(x)51wx52wxdx . (18) Для вариационно-разностной реализации (18) введем на области стержня равномерную сетку юі ={(xi = іЛ), i = 0,1,...,n} на отрезке [0,l], где x = xi -узлы сетки; Л = l / n - шаг сетки по направлению оси координат x . Введем вспомогательную сетку с узлами I: ю^^^Л/2 + іЛ, і = 0,1,..., n -1}. Континуальную область в (18) заменим дискретной, а дифференциальные операторы заменим конечноразностными аналогами: V 51wi+1 - 251 wi +51wi-1 52wi +1 - 252wi + 52wi-1 S = i=2 Л2 Л2 = VN 51wi +1 - 51 wi 52 wi+1 -52 wi Sp 4=1 Л Л (19) где Si и S^ - площадки интегрирования: S^ = Л ; S- = Л - во внутренних узлах i области; Si = Л /2 -в узлах, расположенных на контуре; индексами i = 1 и i = n + 3 обозначены законтурные номера узлов конечно-разностной сетки, а индексами i = 2 и i = n + 2 -номера узлов, расположенных на торцах балки. Тождество (19) содержит произведения вариаций 5jwa и 52wa , для а = 1,2, ..., n + 3, в котором учтем граничные условия (5), принимающие вид: 5jw(0) = 5jw2 = 0, 5jw(l) = 5jwn+2 = 0 , 51w1 = 51w3, 51wn+3 = 51wn+1 ; 52w(0) = 52w2 = 0 , 52w(l) = 52wn+2 = 0, 52w1 =52 w3, 52 wn+3 = 51wn+1 . Применив к выражению (19) циклы а = 3,4,..., n +1, j = 3,4, ..., n +1, k = 3,4, ..., n +1, (20) принимая 51wa = 1, при а = j, 0, при а Ф j, 52 wa = 1, при а = k, 0, при а Ф k, получаем коэффициенты матриц А и В тождества (6). Пусть n = 4, тогда тождество (6), с учетом (19)-(20) представится так: 1 2EJ2 +4EJ3 +EJ4 -2EJ3 -2EJ4 ej4 w3 -2EJ3 -2EJ4 EJ3 +4EJ4 +EJ5 -2EJ4 -2EJ5 w4 EJ4 -2EJ4 -2EJ5 EJ4 +4EJ5 +2EJ6 w5 :л2 Na + Nb - Nb 0 - Nb Nb + Nc - Nc 0 - К Nc + Nd wc (21) Здесь: Na - внутреннее усилие в балке между узлами i = 2 и i = 3; Nb - внутреннее усилие между узлами i = 3 и i = 4 , и т. д. 4 Л 134 Математика, механика, информатика Таким образом, правая и левая части тождества (21) сформированы с использованием вариаций функционала Лагранжа, поэтому матрица жесткости и матрица внутренних сил инерции формируются по единой схеме. В левой части системы (21) матрица жесткости симметричная и содержит суммы параметров жесткости, что позволяет выполнять расчеты стержней переменной жесткости, как и должно, быть при интегральном подходе. В правой части (21) матрица внутренних инерционных сил по сравнению с матрицей правой части в (7) симметричная относительно главной диагонали, нули на главной диагонали не исключаются - последнее не препятствует решению проблемы собственных значений матриц. Продемонстрируем методику расчета на примере действия на стержень (рис. 1) осевого ускорения; для наглядности процедуры вычисления нанесем на область стержня конечно-разностную сетку всего с четырьмя внутренними узлами. Согласно (9) внутренние усилия равны: Na = 3qxl / 4 , Nb = qxl /8, Nc = -qxl /8, Nd = -3qxl / 4 . Пусть EJ = const, тогда из (21) получим EJ " 7 -4 1 " w3 -4 6 -4 ■ w4 1 -4 7 w5 qxl " 7/8 -1/8 0 " w3 -1/8 0 1/8 ■ w4 0 1/8 -7/8 w5 (22) В правой части (22) матрица B содержит нуль на диагонали, и определитель матрицы B равен нулю -обратной матрицы не существует. Покажем поиск решения системы (22). Для этого приведем план решения проблемы собственных значений (22) приемами (6) и (8) ниже. Формируем системы тождеств: [ А]{Щ = s [ B ]{w}, [ B]{w} = s* [ A]{w} . Здесь: s=зі r- = q. =. qx'3 * s = EJ EJ 16EJ EJ 16EJ EJ 3xIX 2 3x1 3 3x 4x3 Запишем операторы в программе Maple [14]: s:=evalf(Eigenvals(A, B, ve)), * s :=evalf(Eigenvals(B, A, ve)), и получим: s := [Float(œ) 5,570484985 - 5,570484986], s* := [0, 1795175826 - 2, 910-11 - 0, 1795175830]. Собственные числа равны: s = -5,57 (минимальное), * s = 0,179 (максимальное). Вычислим критическую нагрузку: 16EJ 16EJ on 1ЛEJ 3x =-^ s = -^(~5,57) = -89,12-, 16EJ qx = ,3 * l s = 89,12 J l3 0,179 l3 Распечатаем ve - матрицы векторов, соответствующих собственным значениям: 0.1428571423 1.000000000 -0.1934378348 1.000000000 0.7020222581 -0.7020222576 0.14285714280.1934378353 -1.000000000 здесь третий столбец матрицы {-0,19; - 0,7; -1} соответствует минимальному собственному числу, 1.000000000 0.14285714370.1934378363 0.7020222554 1.000000000 0.7020222581 0.1934378345 0.1428571431 1.000000000 здесь первый столбец матрицы {1; 0,7; 0,19} соответствует максимальному собственному числу. При перестановке матриц программа Maple изменила знаки собственных значений и собственных векторов. Матрицу [ B ] удобнее представлять безразмерной, выражая произведение ql в качестве параметра, что дает s = qlX2 / EJ = ql3 n2 / EJ . Отсюда q = n2EJs /13. Обозначим через M - массу всего стержня, ax - ускорение, действующее в продольном направлении стержня. Тогда qxl = M ax и M = p Sl, где р - плотность материала, S - площадь поперечного сечения стержня. Теперь по значению найденной критической нагрузки qx вычисляется искомое критическое ускорение ax = qxl /M = n2EJs / (l2pS) = = n2EJ / (l2pSs ). Продолжим поиск критического ускорения стержня (рис. 2), задав: l = 18 м ; E = 106 Па; p = 1 кг/м3; диаметр сечения трубки 1 см; толщина стенки 1 мм. С целью исследования сходимости решения [14] выполним расчеты на различных сетках, начиная с n = 8 (рис. 2, а). Результаты приведем в табл. 1. На сетке n = 64 сходимость практически обеспечивается, что видно из графика сходимости ускорения ax . 4 X 135 Вестник СибГАУ. 2014. N 3(55) Низшую форму потери устойчивости покажем на рис. 2, б. Для конечно-разностных сеток * n = 8, 16, 32, 64, 128 соответствующие значения s , qx и ax приведены в табл. 1. График сходимости ускорения ax приведен на рис. 3. В этой задаче достаточно назначать n = 32 . 3 b 4 c 5 d Л ■*-► -X- Л •*-► -Xe -Xf -Xg h X ■ X 11 6 7 l 10 -► б Рис. 2. Модель стержня, заделанного по обоим торцам: а - конечно-разностная модель, сетка n = 8; б - низшая (первая) собственная форма Таблица 1 n * s 4x ax 8 0,202 6,76 • 106п 12,18 16 0,74 7,38 • 106п 13,28 32 -2,92 -7,53 • 106п 13,55 64 11,61 7,56 • 106п 13,609 128 46,39 7,6 •106п 13,625 Рис. 3. График сходимости ускорения ax Рассмотрим вторую задачу. В барабане центрифуги (рис. 4), вращающемся вокруг вертикальной оси с угловой скоростью ю, расположены стержни, прикрепленные к стенкам барабана, направленные в радиальном направлении. Радиус барабана R2 = 0, 6 м; стержни трубчатые, длиной l = 0, 2 м, радиусом по центру стенки r = 2 10-3 м, толщиной стенки t = 1 • 10-3 м. Радиус Rj = R2 -1 = 0,2 м. Модуль Юнга материала р = 780 кг/м . Требуется определить критические скорости вращения ю барабана центрифуги с позиции устойчивости стержней. Площадь каждого стержня S = 2nrt = 4л-10-6 м2; осевой момент инерции Jx = nr 3t = 8 10-12 п м4; радиус инерции ri = VJx / S = 0,00141 м. Гибкость стержня Л= l/ri = 141,4 - с воззрения сопротивления материалов - стержень гибкий. Выберем начало системы координат Oxy в центре барабана (рис. 4, а). Тогда продольное инерционное внутреннее усилие N(x) в стержне равно N(x) = (pSю2 /2) (R2 - x2), R1 < x < R2. (23) На торце стержня x = R1 известны силовые факторы: Q(R1) = 0, Mx (R1) = 0, N(R1) = 0. На торце x = R1 известны геометрические зависимости - линейные смещения и угол поворота отсутствуют. Приведем матрицу [ B ] к безразмерному виду. Для этого (23) запишем так: N(x) = (r12 - x2 ) / (r^ - R12 ), R1 < x < R2. Теперь, к примеру для n = 5 (рис. 4, б), получим интегральное тождество (6) в виде 1 -2100 -2 5 -4 1 0 1 -4 6 -4 1 <i w3 0 1 -4 6 -4 w4 0 0 1 -4 7 w 0,12 -0,12 0 0 0 -0,12 0,4 -0,28 0 0 0 -0,28 0,76 -0,48 0 0 0 -0,48 1,2 -0,72 0 0 0 -0,72 1,72 w. w2 w 5 . (24) Здесь при безразмерных матрицах A и B имеем собственное число, равное s ^VpS / (r22 -r12 ) . (25) Учитывая Л = l / n = (R2 - R1 ) / n , из (25) получаем искомые значения угловой скорости в зависимости от n : = n j/(R - ^ ( 2R22 - R12 ) s /(p S ) = /(R2 - R1 )yj(2R22 -R12 )/(pSs*). (26) Выполним расчеты на сетках n = 8, 16, 32, 64, 128 . На сетке n = 128 формы потери устойчивости для первых трех низших собственных чисел покажем на рис. 4, в. Первую, вторую и третью угловые скорости материала стержня E = 2 -10 Па, а плотность обозначим, соответственно, символами ю1 (1) „(2) Ч x 1 a а s 136 Математика, механика, информатика (3) , их значения, полученные для различных сеток, приведем в табл. 2. _Таблица 2 n Ю(1) Ю(2) Ю(3) 8 11,64 29,58 45,36 16 11,04 28,56 45,80 32 10,75 27,92 45,16 64 10,62 27,59 44,69 128 10,64 27,41 44,44 6 0 12 3 4 5 / -Х-9-І «-8-і ( 0 X 0 )( 0 a b c d e X R2 - R б Рис. 4. Стержни, вращающиеся в центрифуге: а - барабан, вращающийся вокруг вертикальной оси; б - конечно-разностная модель стержня; в - формы потери устойчивости соответствующие критическим угловым скоростям: <в(1) = 11,64 об/с - сплошная линия; Ю(2) = 27,41 об/с - штриховая линия; <в(3) = 44,44 об/с - штрихпунктирная линия Рис. 5. График сходимости s На графике сходимости (рис. 5) собственного числа s (25) от густоты сетки, видим достаточность сетки n = 64 . Таким образом, в задачах устойчивости стержней приложение вариационно-разностной формулировки краевой задачи исключает особенности, проявляющиеся при действии инерционных нагрузок - несимметричность матриц внутренних сил с вырождением отдельных строк в сравнении с конечно-разностной аппроксимацией дифференциальной формулировки. Примеры показали быструю сходимость решений в зависимости от густоты конечно-разностной сетки. Приложение вариационно-разностной формулировки позволит выполнять расчеты композитных стержней переменной жесткости.

Об авторах

Рашид Альтавович Сабиров

Сибирский государственный аэрокосмический университет имени академика М.Ф. Решетнева

Email: rashidsab@mail.ru
кандидат технических наук, доцент, доцент кафедры технической механики и кафедры летательных аппаратов

Список литературы

Дополнительные файлы