FEA OF SHEAR FLEXIBLE BEAMS

Abstract

The author considers a task of stability in compression 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, where transverse shear strains occur as the basic kinematic nodal variables. The urgency of the model is shown by means of test buckling analysis of beams. Buckling of sandwich beams and thick layered composite under non-classical boundary conditions beams are investigated.

Full Text

При рассмотрении задачи об устойчивости при осевом сжатии балок, податливых при трансверсаль-ном сдвиге, решение выполняется с помощью метода конечных элементов (МКЭ). Математическая модель базируется на теории пластин Рейсснера-Мидлина. Для вывода основополагающих матриц МКЭ используется вариационный подход. Для этого в геометрически нелинейной постановке получено выражение полной потенциальной энергии балочного конечного элемента, узловые компоненты которого включают осредненные по высоте сечения углы трансверсально-го сдвига. Разработанная конечно-элементная модель протестирована на примерах решения классических задач. Ее актуальность показана расчетами на устойчивость монолитных изотропных балок с низкой трансверсальной сдвиговой жесткостью, трехслойных балок с податливым заполнителем, а также балок, изготовленных из однонаправленного композита (ор-тотропная модель). Полученные решения подтверждаются конечно-элементными расчетами, выполненными в среде пакета COSMOS/M в более строгой (двумерной) постановке. Исследуемая конечно-элементная модель допускает описание граничных условий неклассического вида, когда на торцах балки регламентируется поведение углов трансверсального сдвига. Показано, что для некоторых типов граничных условий задача об устойчивости балки адекватным образом может быть решена только при использовании расчетных моделей, подобных описанной в данной работе. Запишем выражение потенциальной энергии деформации в геометрически нелинейной постановке для задачи об изгибе в плоскости балки, податливой при трансверсальном сдвиге [1]: U = 21 I KSX А X + SZ А Z + SXZ AXZ + SZX А ZX )dXdУdZ, (1) 2 (V) где через S обозначены комплексы, связывающие напряжения SX = CTx +TXZ0X , SXZ = TXZ _ CTX®X , SZ ~ ^ Z XZ ®X , SZX ~ ^XZ +^ Z ^ X . 56 2. Личаргин Д. В. Методы и средства порождения семантических конструкций естественно-языкового интерфейса программных систем : дис. ... канд. техн. наук : 05.13.17. Красноярск, 2004. 3. Личаргин Д. В. Порождение дерева состояний на основе порождающих грамматик над деревьями строк // Вестник СибГАУ. 2009. Вып. 4 (25). С. 33-37. 4. Никитин М. В. Лексическое значение слова. М. : Высш. шк., 1983. 5. Селиверстова О. Н. Контрастивная синтаксическая семантика. М. : Наука, 1990. Вестник Сибирского государственного аэрокосмического университета имени академика М. Ф. Решетнева Через А в (1) обозначены деформации нелинейной теории: (2) ex = 4i I ex + -2 ®2 (5) В итоге выражение потенциальной энергии деформации примет вид U = \ I I j{(^11Ax +xxz 0x )Ax + Vxz }dxdydz = (V) = 2II I{A2 +Zxz (exz +0 xA x ))(^dydz. (V) (6) Подставим в (6) выражение для нелинейной продольной деформации из (2): u=2II) f j-4- (ex+2 < I+ + T„ exz +0x I ex + 2 ®2 (7) dxdydz. Заменим здесь линейную продольную деформацию произвольной точки, отстоящей на z от начальной плоскости, ее выражением (3) и перепишем (7) следующим образом: U=21 f f{ Aii (е x+z% x+2 “2 ) + Т x exz + 0x I 8x + z1x + 2“2 (8) dx dy dz. Запишем выражения для внутренних силовых факторов: Nx = II A11 ( 8x + 2 ®2 + z1x I dydz = (S) ( 2 у = B l8 x + "2 ю2 ^+ C1 x, Mx = II A11z I 8x + 2 ®2 + z1 x I dydz = (S) ( 2 У = C l8x + ■J®2 l + D 1 x , (9) (10) Qx = II Gxzexzdydz = KV : (11) (S ) где где ех - линейная часть продольной деформации, вычисляемая по выражению ex = 8x +z ix . (3) Здесь ех - продольная деформация, а хх - изменение кривизны начальной плоскости: Вертикальную деформацию Az учтем в линейной постановке и примем равной нулю. В выражениях для трансверсальных деформаций и напряжений сохраним только линейные члены. Тогда соответствующий комплекс в (1) примет вид Sxz Axz + Szx Azx = Txz (-®x ) + Txz0x = Txzexz . (4) Запишем выражение для продольных напряжений в нелинейной постановке: B = II Andydz ( S ) C = II A11 xdydz ( S ) D = II AJJz2dydz , ( S ) K = II Gxzdydz . ( S ) (12) Интегрируя (8) по площади поперечного сечения с учетом допущений [2] и физических соотношений (9-12), получим: i Nx Л \ 8„ +- + Mx1 x + QxV x + + Qx 0x 2 Л (13) 8„ +- •dx. Следует отметить, что присутствующие в выражении потенциальной энергии деформации нелинейной теории (13) продольное усилие Nx и изгибающий момент Mx (в общем случае) определяются нелинейными физическими соотношениями, а для перерезывающей силы Qx физическое соотношение сохраняет вид, соответствующий линейной теории. Фигурирующее в (13) произведение Qx0x представляет собой проекцию перерезывающей силы на продольную ось балки. Обозначим приведенную продольную силу: N = Nx + Qx (14) Тогда выражение для потенциальной энергии деформации (13), опустив индекс «х», можно переписать в виде выражения: u=21- N + B— + Q 0 2 2 8 + - 2 M + C— ( 2 ( у (15) 1+Qvr dx, где N и M - линейные усилие и момент, определяемые соотношениями N = Be+Ci , M = Ce + Di. Основные кинематические переменные для возмущенного состояния представим в виде us = u0 + u , ws = w0 + w , Vs = V0 + V , (16) где u0, w0, y0 - продольное перемещение, прогиб и угол трансверсального сдвига исходного (докрити-ческого) состояния, u, w, у - соответствующие параметры возмущения. Запишем выражения потенциальной энергии деформации (15) для исходного (докритического) состояний и для возмущенного состояния и вычтем из второго первое: 57 Математика, механика, информатика AU — Us - U0 — 2 N + В— + ( + Q)0 x (17) — +(Mо + M )x+(Qo + Q )ч/г dx. Для задачи об устойчивости балки при осевом сжатии последнее выражение сводится к виду n2 Л_. С ... + Kц>2 Г dx . (18) AU — 1 f I No (f )’ + - Су dx 2 С w cX2 ô= dw dx xv k11 k12 0 k14 k15 0 2 2 k k23 4 2 k 5 2 k 2 k k33 0 k35 k36 4 4 k k45 0 с и м м k55 6 5 k 6 6 k Re = r11 r12 0 r14 r15 0 r22 0 2r 4 25r 0 0 0 0 0 r44 54r 0 с и м м r55 0 (20) k — 12D k — k — k — k — 6D 11 _ ,3 ’ 12 _ 15 _ ft51 _ ft21 _ ,2 l3 l2 12D 4D k14 — k41 — 3 , k22 — i , k23 — k32 — i 6D k24 — k45 — k54 — k42 — ,2 , k25 — k52 — , l2 l D У 2D k26 — k35 — k53 — k62 — i , k33 — 2Kl2 + 6D 6l k36 — k63 — Kl2 -6D k — 12D k — 4D 6l ’ 44 — l3 ’ 55 — l Минимизация функционала (18) позволяет получить разрешающие соотношения для задачи устойчивости. Воспользовавшись МКЭ, получим разрешающие уравнения задачи устойчивости. На первом этапе необходимо сформировать матрицу жесткости элемента балки (для задачи устойчивости). Следует иметь в виду, что фигурирующие в функционале (18) прогиб w и угол сдвига у относятся к возмущению. Рассмотрим прямолинейный двухузловой конечный элемент балки длиной l. Распределение прогибов по длине элемента будем аппроксимировать полиномом третьей степени, а распределение углов транс-версального сдвига - линейной функцией: w(x) — а1 +а2x + а3x2 + а4x3, у (x) — а5 + а6x , где х - координата локальной системы координат элемента, отсчитываемая вдоль продольной оси балки в направлении от первого узла ко второму. Вектор основных кинематических переменных здесь имеет вид D k56 — k65 — l , k66 — 2Kl2 + 6D 6l 6 1 6 1 r" — 5 , r12 — r21 — 12 21 10 , r14 — r41 — - 5 ' - r15 — r51 — 10 2l 1 l 6 r22 — 22 15 2r 4 — 42 — I 10 , 2r — 5r — i 30 : 1 — 4 sr 1 2l T45 — r54 — _ , r55 — ‘ 55 10 15 Рассмотрим балку длиной L, состоящую из n конечных элементов. Для i-го элемента запишем систему уравнений вида (19): Ke® ô.® + No R.® ô.® = 0, где ô.® - вектор узловых неизвестных i-го элемента: Г Vi+1 Г . .<i)= <0 и R <0 Выполняя традиционную для МКЭ процедуру минимизации функционала (18) по компонентам вектора ô, получим однородное матричное выражение: K. ô. + No R. ô. = 0, (19) где матрицы К. и R. имеют следующие структуры и компоненты: К.= K.w и R.<!) - матрицы, аналогичные (20) соответственно. Глобальный вектор узловых неизвестных запишем в виде выражения Л = {ôi Ô2 ... ô,- ... ô„+1}r, (21) где ô, - вектор узловых неизвестных i-го узла: ôi=\w, \ddri -, dx Запишем изменение потенциальной энергии всей балки. Для этого необходимо для каждого элемента получить функционал вида (18), а затем их сложить: AE^—ŸAU, . i—1 (22) Минимизируя (22) по каждой компоненте глобального вектора неизвестных (21), получим следующую систему однородных алгебраических уравнений: К Л +N0Re Л = 0, (23) где КЕ и Re - глобальные матрицы, которые можно получить из соответствующих матриц отдельных элементов К.® и R.®. В математическом смысле система (23) представляет собой обобщенную задачу на собственные значения, решая которую определяют критическое усилие Nup и соответствующую ему форму потери устой- T 58 Вестник Сибирского государственного аэрокосмического университета имени академика М. Ф. Решетнева чивости. Следует иметь ввиду, что в однородной СЛАУ (23) необходимо учесть граничные условия. Выполним тестовое решение для стальной (Е = 210 ГПа, д = 0,3) балки длиной 1 м, с прямоугольным поперечным сечением (h = 1 см, b = 1 мм). Результат сравним с классическим решением, определяемым по формуле Эйлера: P = n2 л2 EJ (24) где n - номер моды, момент инерции сечения вычисляется по формуле J = bh3/12. Сначала в нашем КЭ расчете зададим граничные условия свободного опирания, когда углы сдвига у на торцах не фиксируются: 0 и x = L: w = 0. (25) Номер значения Авторское решение Аналитическое решение Решение в пакете COSMOS/M 1 172,690 172,718 172,667 2 690,429 690,872 690,058 3 1552,218 1554,462 1550,35 упругой линии балки. Таким образом, необходимо сделать следующий вывод, имеющий практическое значение при решении задачи устойчивости балок с низкой трансверсальной жесткостью по моделям, учитывающим угол трансверсального сдвига: при выборе величины критического усилия из списка собственных значений задачи необходимо анализировать соответствующие им формы потери устойчивости, получаемые, как известно, на основе собственных векторов. Таблица 2 Собственные значения задачи об устойчивости шарнирно-опертой балки с уменьшенной трансверсальной сдвиговой жесткостью при x В результате решения системы (23) для 50-узловой модели с граничными условиями (25) получено критическое усилие сжатия N^, равное 172,69 (Н). Критическая нагрузка, вычисленная для балки тестового примера по формуле (24), оказалась равной 172,72 (Н). Почти идеальное совпадение численного решения с аналитическим подтверждает правильность рассматриваемого здесь алгоритма. Можно сделать более строгий анализ полученных нами результатов. Для этого сравним величины нескольких первых собственных значений нашей задачи устойчивости с аналогичными величинами, полученными по формуле Эйлера для старших мод. Результаты нашего решения и значения, вычисленные по формуле (24), приведены в табл. 1, где также представлены данные, полученные при решении задачи устойчивости в пакете COSMOS/M, в расчетной модели которой использовались элементы PLANE2D. Таблица 1 Собственные значения задачи об устойчивости шарнирно-опертой балки Номер Собственные значения 1 33,120 2 131,507 3 169,99 4 649,169 5 1357,860 Выполним решение задачи устойчивости балки при осевом сжатии, уменьшив в 100 раз (по сравнению с предыдущим расчетом) жесткость балки при трансверсальном сдвиге. Все остальные параметры жесткости, размеры и граничные условия оставим без изменения. Результаты решения по нашей модели представлены списком первых пяти собственных значений задачи (табл. 2). Третье, четвертое и пятое значения (табл. 2) на самом деле соответствуют первым трем формам потери устойчивости балки. Первые два значения (33,120 и 131,507) соответствуют изменению формы за счет изменения углов наклона сечений (очень малых абсолютных величин) при сохранении прямолинейности Рассмотрим трехслойную балку с податливым заполнителем. Решим задачу об ее устойчивости при осевом сжатии. Назначим следующие геометрические параметры: длина - 1 м, толщина несущих слоев - 1 мм, толщина слоя заполнителя - 3 см. Материал несущих слоев - сталь (Е = 210 ГПа, д = 0,3). Изотропный материал заполнителя обладает модулем сдвига Gзап = 68 МПа. Результаты решения задачи на собственные значения (23), (25) при ширине сечения в 1 мм приведены в табл. 3. Здесь, как и в расчете монолитной балки, собственные значения, соответствующие формам общей потери устойчивости, определены при анализе собственных векторов системы. При этом критическое усилие сжатия оказалось равным 728,8 (Н). Рассмотрим балку с различными вариантами защемления торцов. Сначала выполним тестовое решение для изотропной (сталь) монолитной балки. Критическое усилие оказалось равным 690,429 (Н). Этот результат почти совпадает с аналитическим решением, вычисленным по формуле P^ = 4 п2 EJ/L2. Таблица 3 Собственные значения задачи об устойчивости трехслойной шарнирно-опертой балки Номер Собственные значения 1 54,912 2 209,67 3 438,433 4 708,916 5 728,802 (1 основное) 6 991,453 7 1263,991 8 1512,876 9 1608,109 (2 основное) 10 1731,201 11 1916,448 12 2053,940 (3 основное) 59 Математика, механика, информатика Полученные результаты хорошо согласуются с данными расчета, выполненного в пакете COSMOS/M (табл. 4) с помощью плоской модели, в которой задействованы элементы PLANE2D. При этом обнаружено качественное совпадение по формам. Некоторое расхождение по собственным числам объясняется различием КЭ моделей, связанным прежде всего с мерностью задачи: напомним, что решение в пакете COSMOS/M выполнено с помощью двумерной модели. Таблица 4 Собственные значения задачи об устойчивости изотропной балки, защемленной по обоим торцам Номер Авторское Решение в пакете значения решение COSMOS/M 1 690,43 673,36 2 1411,31 1376,87 3 2756,39 2684,58 Еще одно тестовое решение выполнено для стальной жестко защемленной балки, представляющей собой узкую пластину: ширина сечения - 50 мм, высота - 10 мм. В этом случае результаты можно сравнить с расчетом, выполненном в пакете COSMOS/M, в котором пространственная конструкция моделируется элементами слоистой оболочки SHELL4L. При этом граничные условия более адекватны тем, что используются в нашей модели. Полученные собственные значения двух решений (табл. 5) по всем трем формам имеют расхождение менее 1 %, что подтверждает точность нашей модели. Решена модельная задача об устойчивости жестко защемленной по обоим торцам балки с уменьшенной в 200 раз (по сравнению с изотропной) трансверсаль-ной жесткостью. При этом критическая нагрузка уменьшилась до 612,214 (Н). Собственные значения, соответствующие второй и третьей формам, оказались равными 1096,364 и 1822,815 (Н) соответственно. Таблица 5 Собственные значения задачи об устойчивости полосы, защемленной по обоим торцам Номер значения Авторское решение Решение в пакете COSMOS/M Расхождение, % 1 34521,44 34763 0,69 2 70565,55 71093,4 0,74 3 137819,65 139083 0,908 Результаты решения задачи на устойчивость трехслойной балки длиной 1 м, со стальными несущими слоями толщиной по 1 мм, с податливым заполнителем (Езап = 2 ГПа) толщиной 30 мм приведены в табл. 6 и 7. Рассматривались два вида граничных условий: строгое защемление (у = 0, табл. 6) и защемление со свободным сдвигом (у ф 0, табл. 7). Параллельно задача решалась в пакете COSMOS/M, где трехслойная конструкция моделировалась элементами SHELL4L. Выявлено хорошее совпадение результатов по нашей модели и решением в пакете COSMOS/M для случая строгого соответствия граничных условий. Такое соответствие можно обеспечить только для полного защемления (табл. 6). Второй вид защемления (со свободным сдвигом) в случае весьма податливого заполнителя можно моделировать в пакете COSMOS/M условиями шарнирного опирания. Данные табл. 7 свидетельствуют о возможности такой «подмены». Таблица 6 Собственные значения задачи об устойчивости трехслойной балки с полным защемлением торцов Номер значения Авторское решение Решение в пакете COSMOS/M Расхождение, % 1 173757 177339 2,02 2 306423 304751 0,55 3 501557 486576 3,08 Таблица 7 Собственные значения задачи об устойчивости трехслойной балки, на защемленных торцах которой допускается сдвиг Номер значения Авторское решение Решение в пакете COSMOS/M 1 52523,9 49892,1 2 176399 176796 3 336925 334228 Рассмотрим задачу об устойчивости балки, изготовленной из однонаправленного композита, например, углепластика, у которого следующие механические свойства: Е = 180 ГПа, G = 5 ГПа. Выполним серию расчетов для различных вариантов граничных условий. Сначала торцы балки будем полагать жестко защемленными. Результаты в виде первых трех собственных значений приведены в табл. 8. Таблица 8 Собственные значения задачи об устойчивости композитной балки с полным защемлением торцов Номер Авторское Решение в пакете значения решение COSMOS/M 1 586,95 556,25 2 1187,66 1121,5 3 2287,01 2143,23 Здесь же для сравнения представлены результаты решения, выполненного в пакете COSMOS/M, где в двумерной модели балки использовались элементы PLANE2D, допускающие задание свойств ортотроп-ного материала. В результате аналогичного решения, выполненного для шарнирно-опертой балки со свободным сдвигом на торцах (табл. 9), обнаружено хорошее совпадение собственных значений, определенных по нашей модели, с числами, полученными при КЭ расчете в пакете COSMOS/M, где рассматривалась двумерная модель с элементами PLANE2D. 60 Вестник Сибирского государственного аэрокосмического университета имени академика М. Ф. Решетнева Таблица 9 Собственные значения задачи об устойчивости шарнирно-опертой композитной балки Номер моды Тестируемое решение Решение в пакете COSMOS/M 1 147,716 147,538 2 586,955 584,162 3 1306,212 1292,53 При расчете балки из однонаправленного углепластика со вторым вариантом защемления, когда угол сдвига на торцах не фиксируется, определены собственные значения, первые два из которых оказались равными 71,208 и 321,916. Соответствующие им формы показаны на рисунке. Обратим внимание на то обстоятельство, что ни при шарнирном опирании торцов, ни при их полном защемлении среди форм общей потери устойчивости не было одновременно двух форм с одной полуволной. Критическая форма (рисунок, а) по своему виду соответствует форме потери устойчивости шарнирно-опертой балки. Такая картина возможна, как мы отмечали ранее, при весьма низкой трансверсальной жесткости балки. На графиках распределения прогибов второй формы (рисунок, б) явно заметно влияние защемления торцов. В заключение сформулируем выводы: 1. Получен энергетический функционал для решения задачи об устойчивости балок с низкой трансвер-сальной жесткостью. 2. Разработана конечно-элементная модель податливой при сдвиге балки, вектор узловых неизвестных которой включает углы трансверсального сдвига. Сформированы матрица жесткости и матрица геометрической жесткости для соответствующего балочного элемента. 3. В результате проведенного численного исследования на примере изотропных, ортотропных композитных и трехслойных балок показана актуальность разработанной модели при проведении расчета на устойчивость конструкций с низкой трансверсальной сдвиговой жесткостью, а также при учете граничных условий неклассического вида