Modal finite element analysis of shear flexible beams

Abstract

The author considers a problem of free oscillations of a beam with low value of transverse shear stiffness. On the base of Reissner-Midlin plate theory, the finite element of a shear flexible beam is developed. There are transverse shear strains among the basic kinematic nodal variables. The urgency of the model is shown by test modal analysis of beams. Vibrations of sandwich beams and thick layered composite beams under non-classical boundary conditions are investigated.

Full Text

Композиционные материалы часто используются в производстве ракетно-космической техники, поскольку они, обладая высокими удельными механическими свойствами, позволяют изготавливать конструкции с большой степенью весового совершенства. Для композитов характерна низкая трансверсальная сдвиговая жесткость. Эту особенность поведения композитов необходимо учитывать при проектировании и расчете. Математическая модель при этом усложнится по сравнению с той, что используется при расчете конструкций из традиционных материалов. Рассматривается задача о собственных колебаниях балок, податливых при трансверсальном сдвиге. Решение выполняется с помощью метода конечных элементов (МКЭ). Математическая модель базируется на теории пластин Рейсснера-Мидлина [1; 2]. Для вывода матриц жесткости и инерции используется вариационный подход. Для этого выражение полной энергии конечного элемента колеблющейся балки варьируется по узловым компонентам, среди которых присутствуют осредненные по высоте сечения углы трансверсального сдвига. Разработанная конечноэлементная модель протестирована на примерах решения классических задач. Ее актуальность показана расчетами, выполненными для монолитных изотропных балок с низкой трансверсальной сдвиговой жесткостью, для трехслойных балок с податливым заполнителем, а также для балок, изготовленных из однонаправленного композита (ортотропная модель). Частоты и формы собственных колебаний подтверждаются конечно-элементными решениями в среде пакета COSMOS/M, где эти задачи рассматривались в более строгой (двумерной) постановке. Исследуемая конечно-элементная модель допускает описание граничных условий неклассического вида, когда на торцах балки регламентируется поведение углов трансверсального сдвига. Показано, что для некоторых типов граничных условий задача о собственных колебаниях балки не может быть решена при использовании классических расчетных моделей. Получим разрешающие уравнения для задачи о свободных поперечных колебаниях балки с низкой трансверсальной сдвиговой жесткостью. Задачу будем решать с помощью МКЭ, поэтому сначала необходимо записать выражение энергетического функционала. Полная энергия колеблющейся балки складывается из потенциальной энергии деформации и кинетической энергии движения: E = Ü + T . (1) Здесь + Txzêxz ) dxdydz .. Ü = 21 J J( +^z (V) T= 1J J Ip(x, y, z) v2 dxdydz , (2) 2 (V) где p - плотность материала в точке с координатами (x, y, z); v - скорость движения этой точки тела; V - объем тела. Последовательно применяя гипотезы о поведении балки с низкой трансверсальной сдвиговой жесткостью [3], вместо (1) получим 1 Г f -кт —u , ,d 9 ^ ^JI N-т + M— + Qy I dx- 20 V dx dx J 4 J k(u2+w )+2cpu fy-î J+°p(y-—xJ }x x dx = 0 . (3) где N, M, Q - внутренние балочные усилия (N - продольная сила, М - изгибающий момент, Q - перерезывающая сила), вычисляемые по формулам N = B,—uL + C dx f ,1 J2 'l d y d w dx dx2 M = C^u + D dx f J j2 'l d y d w dx dx2 Q = K y, (4) где Bp, Ср, Dp - параметры, характеризующие инерционные свойства b 2 h-s b_ 2 h-s Bp = J J pdydz, Cp = J J pzdydz, _b -s 2 b_ 2 h-s Dp = J J p z2 dydz. (5) _ b -s 2 69 Вестник СибГАУ. № 1(47). 2013 Здесь h и b - высота и ширина сечения балки соответственно. Перейдем теперь непосредственно к формированию матрицы жесткости балочного конечного элемента. Будем полагать, что прогиб изменяется по кубическому закону вдоль длины элемента, а угол сдвига у и продольное перемещение и - по линейному закону: w (x) = a1 +а2x + а3x2 +а4x3, у (x) = а5 + а6x , u (x) = а7 +а8x. (6) Введем в рассмотрение вектор основных кинематических переменных я J dw , о = \ w — у u f . dx dx dx 26 ^36 77 D D D D =12 F , k12 = 6 “2 , l k15 =-1273, k16 = 612 D D k C D = 4— l , k23 =--, 23 l 24 = l ’ k25 = -6 if D D CD Kl 1 C = 2— l , k27 = "y > k28 = —, k33 = —+ ■ l 33 l , k34 — 3 34 l = D D Kl C B = l ’ koT — + 37 ! 6 , k38 = i л % = i C C k B D =7 ’ k47 =- 7 ’ =- 7 ’ k55 = 12 ~3 D D D C = -6 ? ^ k66 =4 7, k67 = - 7 ’ k68 =- "7 = D = l 13BP Kl V C k 3 , k78 = I •> k88 l2 + 42DP = B = i ’ 11BP l2 Dp Dp p p 1 о • _ P , P _ p 35l 210 10 Cp p 3(3-ßp 12 - 28^pj _ 13Bp i2 D S14 = 2 , S15 = , = +" 43 2 Dp 2 70l Dp S = C 2 , 18 2 46 420 10 s17 = — , s17 =—p , s18 = —- , s22 = —^Bp ! +—Dp l 17 „ > 17 „ . 18 „ > 22 105 p 15 p S23 = S24 = Dpi 12 1 3 1 s26 =--Bp l3--Dp l 26 140 p 30 p Cpl 12 S27 = S25 = Dp l S33 = S37 = Dp l S34 = Cp 1 12 Dp 13Bp l2 420 S28 = S35 = " S36 = (7) Вектор узловых кинематических параметров произвольного элемента балки представим в следующем виде: ôe = Jw1 f У1 U1 W2 fУ 2 U2} . (8) 66 S77 = .Dp! 6 ’ Cp l S38 = — , Bp l s44 = —— 44 3 , S45 Cp i Cp i Bp i 13Bp l2 —— , s47 12 47 = , , S48 6 = 6 , S55 =3 11Bp l2 - Dp <N 1 = 7 5 S 210 10 , S58 — bp i3 105 p +—?P l , 15 p = DPl , S68 = S67 12 _Dp 10 Cp! 12 Dp l Cp Cp Cl 12 Dp! 3 S78 = Cp! 3 S88 = Bp! 3 Следуя вариационному подходу к МКЭ, получим систему уравнений для произвольного балочного элемента, которая в матричном представлении имеет следующий вид: Ke ôe - œ2 Se ôe = 0, (9) где К и Se - матрица жесткости и матрица инерци-альных параметров балочного элемента балочного элемента со следующими ненулевыми компонентами: С помощью матриц Ke и Se формируются глобальные матрицы жесткости и инерции Ks и Ss, фигурирующие в глобальных уравнениях состояния К Л - œ2 Ss Л = 0, (10) где глобальный вектор узловых неизвестных имеет вид Л = {Ô1 02 ... Ôi ... 0Я+1>Т, ô; - вектор узловых неизвестных i-го узла (11) dw dx У г иг Описанный алгоритм вывода матрицы жесткости и матрицы инерции реализован в специальной программе формирования фундаментальных матриц и векторов теории МКЭ, предназначенных для расчета балок, податливых при трансверсальном сдвиге [4]. Выполним тестовое решение для стальной (Е = 210 ГПа, д = 0,3, р = 7 800 кг/м3) балки длиной 1 м, с прямоугольным поперечным сечением (h = 1 см, b = 1 мм). Будем полагать, что координатная ось X проходит через центры тяжести сечений. В этом случае смешанные жесткости (С), фигурирующие в матрице жесткости элемента К ,и параметр инерции (Ср), присутствующий в матрице Se, равны нулю. Сначала в нашем КЭ расчете зададим граничные условия свободного опирания, причем углы сдвига у на торцах не фиксируются: при x = 0: w = 0, и = 0, при x = L: w = 0. (12) Результат сравним с классическим решением, определяемым по формуле (13) где J - момент инерции сечения. С помощью программы [5] решена обобщенная задача на собственные значения для 50-узловой моде- 3 3 2 2 2 S 46 56 2 Т 22 46 12 2 70 Математика, механика, информатика ли стальной балки длиной 1 м с прямоугольным поперечным сечением (h = 1 см, b = 1 мм) и получено значение первой круговой частоты, равное 147,811 (рад/с). Вычисленное по формуле (13) значение равно 147,833 (рад/с). Форма колебаний при этой частоте показана на рисунке. Представлены первые пять частот собственных колебаний стальной балки, определенные по нашей модели, вычисленные по формуле (13) и найденные в результате решения в пакете COSMOS/M с использованием элементов PLANE2D (табл. 1). Почти идеальное совпадение по частотам подтверждает правильность рассматриваемого здесь подхода. Совпадение имеет место и по формам колебаний. Отметим, что результаты модального анализа в пакете COSMOS/M выдаются в виде циклической частоты f, измеряемой в Гц (с-1) и связанной с круговой частотой формулой œ = 2 п f Первая форма колебаний балки Таблица 1 Частоты собственных колебаний для первых пяти мод Важно отметить, что в данном случае частоты первых мод колебаний являются первыми в списке собственных значений задачи, т. е. наименьшими. Выполним модальный анализ, уменьшив в 100 раз (по сравнению с предыдущим расчетом) жесткость балки при трансверсальном сдвиге. Все остальные параметры жесткости, размеры и граничные условия оставим без изменения. Результаты решения по нашей модели представлены списком первых пяти собственных значений задачи, расположенных во втором столбце (табл. 2). Здесь же для сравнения приводится решение (табл. 2, третий столбец), выполненное в пакете COSMOS/M, в котором для учета уменьшенного параметра трансверсальной жесткости задействуется опция ортотропных свойств материала (для элементов PLANE2D). При этом фактически уменьшается соответствующий модуль сдвига в списке механических свойств материала. Вновь отмечаем хорошее совпадение результатов. Подчеркнем необходимость учета трансверсального сдвига в модальном анализе. Сравнение соответствующих данных (см. табл. 1 и 2) показывает заметное расхождение частот собственных колебаний, определенных с учетом трансверсального сдвига, и теми же величинами, вычисленными по классической формуле, где этот эффект не принимается во внимание. Выполнено решение задачи о собственных колебаниях балки с шарнирно-опертыми концами, на которых зафиксирован угол трансверсального сдвига. В этом случае рассматривается разрешающая система алгебраических уравнений со следующими граничными условиями: при x = 0: w = 0, u = 0, у = 0, при x = L: w = 0, у = 0. (14) Результаты решения (первые пять частот) приведены в четвертом столбце (табл. 2). Они показывают, что для рассматриваемой здесь шарнирно-опертой балки (с данным набором геометрических и механических параметров) расхождение в собственных значениях, определенных по двум моделям (с фиксацией угла у на торцах и без такового), невелико. Тем не менее качественное расхождение - верное: фиксация угла у на торцах делает балку более жесткой, что отражается в повышении величин собственных частот. Таблица 2 Частоты колебания шарнирно-опертой балки Номер моды Торцы со свободным сдвигом Решение в пакете COSMOS/M Торцы с фиксированным сдвигом 1 146,271 146,058 146,308 2 567,477 564,397 568,005 3 1 218,091 1 204,820 1 220,416 4 2 040,959 2 006,992 2 047,117 5 2 982,133 2 915,593 2 994,430 Параметры инерции балки в данном случае для прямоугольного сечения с шириной b и высотой h определяются следующим образом: b h3 Bp=pbh , Cp= ^ Dp=p —. Рассмотрим трехслойную балку с податливым заполнителем. Решим задачу о ее собственных колебаниях. Назначим следующие геометрические параметры: длина - 1 м, толщина несущих слоев - 1 мм, толщина слоя заполнителя - 3 см. Материал несущих слоев - сталь (Е = 210 ГПа, д = 0,28, р = 7 700 кг/м3). Изотропный материал заполнителя обладает модулем сдвига Gsan = 2,9 МПа. С помощью программы [6] получены результаты решения задачи на собственные Номер моды Тестируемое решение Аналитическое решение Решение в пакете COSMOS/M 1 147,811 147,833 147,809 2 590,982 591,332 590,956 3 1 328,730 1 330,497 1 328,605 4 2 359,758 2 365,328 2 359,361 5 3 682,269 3 695,826 3 681,293 71 Вестник СибГАУ. № 1(47). 2013 значения при ширине сечения в 1 мм (табл. 3). Параметры инерции трехслойной балки вычисляются по формулам Bp =(н hn +Pз h )b , Cp = ^ Dp = (p„ (3 - hl ) + p з h3 ), где рн и рз - плотности материалов несущих слоев и заполнителя; hK и hs - толщины несущих слоев и слоя заполнителя; H - общая толщина пакета, H = 2 hK + h3. Таблица 3 Частоты колебания трехслойной шарнирно-опертой балки Для сравнения приводится результат (табл. 3, третий столбец) решения задачи в пакете COSMOS/M, где трехслойная конструкция моделировалась элементами PLANE2D (заполнитель) и BEAM2D (несущие слои). Хорошее совпадение частот собственных колебаний (табл. 3), вычисленных с помощью двух разных моделей, подтверждает правильность разрабатываемого нами подхода. Здесь же приводятся результаты расчета для модели с фиксированным сдвигом на торцах (см. табл. 3, четвертый столбец). Небольшое отличие частот для балок с фиксацией сдвига и без такового характерно для граничных условий шарнирного опирания. Выполнен модальный расчет трехслойной конструкции в пакете COSMOS/M с моделированием балки элементами слоистой пластины SHELL4L. На торцах -условия шарнирного опирания: слева - неподвижный шарнир, справа - подвижный в направлении продольной оси. Ширина сечения балки - 50 мм. Остальные размеры и механические свойства материалов такие же, как в балке предыдущего расчета. Результаты расчетов (табл. 4) сравниваются с данными, определенными по нашей модели. Как видим, имеется заметное расхождение, что свидетельствует о том, что в модели пакета COSMOS/M не учитываются рассматриваемые нами эффекты. В данном случае большего доверия заслуживают наши результаты, поскольку они почти точно совпадают с решением, полученным в том же пакете COSMOS/M при строгом моделировании поведения заполнителя (см. табл. 3). Рассмотрим трехслойную балку с граничными условиями жесткого защемления торцов. при x = 0 и при x = L: w = 0, u = 0, у = 0, dw/dx = 0. (15) Сохраним геометрические параметры предыдущей задачи: длина - 1 м, толщина несущих слоев - 1 мм, толщина слоя заполнителя - 3 см, ширина сечения - 1 мм. Материал несущих слоев - сталь (Е = 210 ГПа, д = 0,28, р = 7 700 кг/м3). Модуль сдвига изотропного материала заполнителя примем Gs^ = 290 МПа. Таблица 4 Частоты колебания трехслойной полосы Номер моды Торцы со свободным сдвигом Решение в пакете COSMOS/M (SHELL4L) 1 20,964 22,227 2 43,292 46,028 3 65,338 69,487 4 87,306 92,829 5 109,243 116,097 Рассмотрим результаты решения задачи (табл. 5). Для сравнения приводятся частоты (см. табл. 5, третий столбец), вычисленные с помощью пакета COSMOS/M, где трехслойная конструкция моделировалась элементами PLANE2D (заполнитель) и BEAM2D (несущие слои). Защемление торцов обеспечено фиксацией угловых узлов на краях. Вновь следует отметить весьма хорошее совпадение и в количественном (табл. 5) и в качественном отношении результатов, полученных по нашей модели и в пакете COSMOS/M. Выполнен модальный расчет для балки с защемленными торцами, на которых, однако, освобожден угол сдвига. Значения частот собственных колебаний показаны в предпоследнем столбце табл. 5. Они существенно расходятся с соответствующими величинами, вычисленными для жестко защемленной балки (табл. 5, второй столбец). Следовательно, в модальном расчете балки со вторым видом защемления для точного определения частот и форм собственных колебаний требуется рассматриваемая в нашей работе модель (постановка). Для сравнения выполнено решение задачи для балки с шарнирным опиранием торцов, где, кроме того, допущен сдвиг (см. табл. 5, последний столбец). Как видим, результаты отличаются от тех, что получены для защемленной балки со свободным сдвигом на торцах. Последнее решение задачи об определении частот собственных колебаний защемленной (со свободным сдвигом на торцах) трехслойной балки выполнено для весьма податливого заполнителя Gзап = 2,9 МПа. Получены следующие величины пяти первых частот: f = 21,262 Гц, f2 = 43,833 Гц, f = 66,136 Гц, f4 = 88,365 Гц, f5 = 110,567 Гц. Сравнивая эти значения с частотами колебания шарнирно-опертой балки со свободным сдвигом на торцах (см. табл. 4), можно заключить, что в случае достаточно малой жесткости балки при трансверсальном сдвиге защемление торцов со свободным сдвигом можно моделировать шарнирным опиранием (тоже со свободным сдвигом). В заключение выполним модальный анализ монолитной балки, изготовленной из однонаправленного композита, например углепластика, у которого следующие актуальные для данной задачи механические Номер моды Торцы со свободным сдвигом Решение в пакете COSMOS/M Торцы с фиксированным сдвигом 1 20,99 20,388 21,220 2 43,304 42,278 43,805 3 65,346 64,473 66,119 4 87,312 87,22 88,341 5 109,248 111,049 110,555 72 Математика, механика, информатика свойства: Е = 180 ГПа, G = 5 ГПа и плотность таты этих вычислений в виде первых пяти частот собр = 1500 кг/м3. Выполним серию расчетов для различ- ственных колебаний (табл. 6). ных вариантов граничных условий. Получены резуль- Таблица 5 Частоты колебания трехслойной балки, защемленной по торцам Номер моды Жесткое защемление торцов Решение в пакете COSMOS/M Защемление со свободным сдвигом на торцах Шарнирно-опертые торцы со свободным сдвигом 1 134,105 130,53 77,563 70,122 2 309,411 298,818 254,025 245,059 3 518,337 497,572 477,927 467,133 4 742,287 710,749 716,150 703,270 5 972,990 931,689 956,270 941,120 Таблица 6 Частоты колебания многослойной балки при различных вариантах закрепления торцов Номер моды Шарнирно-опертые торцы со свободным сдвигом Шарнирно-опертые торцы с фиксированным сдвигом Защемление торцов Авторская модель COSMOS/M фиксированный сдвиг свободный сдвиг 1 49,598 49,587 49,599 111,787 96,267 2 197,494 197,335 197,521 305,515 269,699 3 441,062 440,279 441,197 592,381 533,552 4 776,122 773,732 776,538 966,163 884,950 5 1 197,208 1 191,62 1 198,192 1 420,876 1 319,846 Как и в предыдущих примерах, здесь рассматривались следующие варианты граничных условий: классическое шарнирное опирание со свободным транс-версальным сдвигом на торцах, шарнирное опирание с фиксацией сдвига на торцах, жесткое защемление торцов, защемление со свободным сдвигом на торцах. Здесь же приведены данные расчета, выполненного для балки с классическим шарнирным опиранием в пакете COSMOS/M, где монолитная композитная конструкция, колебания которой происходят в одной плоскости, моделируется элементами PLANE2D. Как следует из представленных значений (см. табл. 6), имеется хорошее совпадение данных нашего расчета с результатами вычислений, выполненных в пакете. Напомним, что наше решение соответствует одномерной балочной модели, а расчеты в пакете COSMOS/M проводились с помощью двумерной модели, чтобы адекватно учесть поведение конструкции с низкой трансверсальной сдвиговой жесткостью. Анализ представленных данных (см. табл. 6) подтверждает ранее сделанный вывод о том, что в общем случае модальный расчет композитных балок с неклассическими граничными условиями (шарнирно-опертые торцы с фиксированным сдвигом и защемление со свободным сдвигом на торцах) приводит к результатам, отличным от тех, которые получаются при классических вариантах этих граничных условий. Это обстоятельство подчеркивает актуальность разработанной модели. В заключение сформулируем выводы. 1. Получен энергетический функционал для решения задачи о собственных колебаниях балок с низкой трансверсальной жесткостью. 2. Разработана конечно-элементная модель податливой при сдвиге балки, вектор узловых неизвестных которой включает углы трансверсального сдвига. Сформированы матрица жесткости и матрица инерции для соответствующего балочного элемента. 3. В результате проведенного численного исследования на примере расчета изотропных, ортотроп-ных композитных и трехслойных балок показана актуальность разработанной модели при проведении модального анализа конструкций с низкой трансвер-сальной сдвиговой жесткостью, а также при учете граничных условий неклассического вида.