CC BY
185
А. Х. Валиуллин
БАЛОЧНЫЕ ЭЛЕМЕНТЫ ВЫСОКОГО ПОРЯДКА
ДЛЯ РАСЧЕТА РАМ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ
Ключевые слова: конечный элемент, локальная матрица жесткости, глобальная матрица
жесткости, расчет сложных рам.
На основе вариационного принципа Лагранжа получены конечные элементы в форме балок с аппроксимацией прогиба полиномами четвертой и пятой степени. Получены матрицы жесткости элемента. Показан порядок формирования глобальной матрицы жесткости и вектора правой части разрешающей системы уравнений. Предлагаемые элементы позволяют выполнить расчет НДС сложных статически неопределимых рам, нагруженных как сосредоточенными, так и распределенными силами постоянной и линейнопеременной интенсивности.
Key words: finite elements, local stiffnes matrix, global stiffnes matrix, galculation of the
complex frames.
Finite elements in the form of beams with approximation of normal displacements by polynoms of fourth or fifth degrees are received on the basis of Lagrange variation principle. Stiffness matrixes of these elements are obtained. The procedure of forming of the governing equations is shown. Offered elements allow to execute calculation of the complex frames loaded both concentrated and allocated external forces.
При расчете стержневых систем обычно используется элемент третьего порядка [1], когда прогиб аппроксимируется полиномом третьей степени:
w = Wj(l -3^2 + 2Ч1)+ /,9Д1 -Ч+ + Wk (зЧ2 -2Ч3)+ ¡АЧ1 (-1). (1)
Здесь обозначено (рис. 1): ¡, - длина i - го элемента, Ч = x/lj, x, z - локальные
координаты, Uj, Wj, j, uk, wk, Qk - локальные перемещения узлов элемента.
Рис. 1 - Балочный двухточечный конечный элемент третьего порядка
Элемент обеспечивает достаточную точность и широко используется, однако в случаях немонотонного изменения прогиба в пределах одного элемента обнаруживается погрешность. Так, например, при представлении одним элементом балки, свободно опертой по концам и несущей равномерно распределенную нагрузку, ее прогиб посредине получается равным ц!4 /(96ЕЛ), тогда как точное решение дает 5ц!4 /(384ЕЛ) = = ц!4 /(76,8EJ). Разница, как видно, значительна; кроме того, в местах начала и конца распределенной нагрузки наблюдаются невязки в эпюрах изгибающего момента величиной порядка ц!2 /12. Правда, эти недостатки можно устранить измельчением и увеличением числа элементов. Если ту же балку разбить на два элемента, результат практически совпадает с точным решением. Причина наблюдаемых отклонений состоит в недостаточности степени полинома, принятого для аппроксимации прогиба.
Здесь предлагаются элементы четвертого и пятого порядка, которые можно использовать при расчете сложных рамных и балочных конструкций, нагруженных как сосредоточенными, так и распределенными по длине силами. При этом нет необходимости в дроблении элементов, в качестве последних можно принимать целиком участки конструкции, границами которых могут быть места приложения сил, углы рамы, опоры, промежуточные шарниры и т. п.
Получим аппроксимационную формулу, а затем и матрицу жесткости элемента для полинома четвертой степени, а для полинома пятой степени приведем только результаты, полученные аналогично.
Рис. 2 - Балочный трехточечный конечный элемент
Элемент (рис.2) длиной !,■ имеет три узловых точки у, т, к; начало местных координат принимаем в середине элемента:
- !,/2 < х < !,/2 (- 0,5 < £ < 0,5).
Представим прогиб в сечении х полиномом четвертой степени:
w = а1х4 + а2 х3 + а3 х2 + а4х + а5, (2)
угол поворота
dw 3 2
д =— = 4а,х3 + 3а2 х2 + 2а3 х + а4. (3)
dх 1234 V 7
Коэффициенты а найдем из следуюших граничных условий: при х = 0 w = wm, при х = -0,5!, w = Wj, д = ду, при х = 0,5!, w = wk, д = дк . Подставив найденные коэффициенты в формулы (2) и (3), получим в окончательном виде выражения прогиба и угла поворота:
w = wj (- 8£4 + 2£3 + 4£2 -1,5£)+ !1 ду (- 2£4 + £3 + 0,5£2 - 0,25^) +
+ Wm(16£4 - 8£2 +1)+ wк(- 8£4 - 2£3 + 4£2 +1,5£)+ !кдк(2£4 + £3 - 0,5£2 - 0,25^); (4)
3 = - [Wj (- 32^3 + 6^2 + 8£-1,5) + ) 3 j (- 8^3 + 3^2 + £- 0,25)+
¡і
+ wm (б4^3 -16^)+ wk (- 32^3 - 6^2 + 8^ +1,5)+¡i3k (8^3 + 3^2 -£- 0,25); . (5)
Для составления матрицы жесткости элемента воспользуемся вариационным принципом Лагранжа [2]. Составим условие экстремума дополнительной энергии:
5Э = 5(W - U ) = 0, (6)
0,51
где SW = JqSwdx +^Fk5wk +^ Mn53n - вариация работы внешних сил: распределен-
- 0,5l k n
ной нагрузки q, сосредоточенных сил Fk и сосредоточенных моментов Mn, sm V Гсл dud _ .d 2wd 2w S V
oU = J EA------------------ou + EJ—2-гSw dx - вариация потенциальной энергии
- 0,5l
dx dx dx 2 dx
деформации (без учета сдвига).
Примем во внимание, что продольное перемещение аппроксимируется линейной функцией
и = и у (0,5 -£) + ик (0,5 + £). (7)
Подставим выражения перемещений и и w и их производных в уравнение (6), предварительно записав их в матричном виде:
^ = 1 (- иу + их )=к 1(4
где [рм] = -1 [-1 0 0 0 0 0 1 0 0], (и}- матрица-столбец перемещений узловых
точек элемента, запишем его в виде строки:
(и}Г = [ и/у 3у ит wm дт ик Wk дк]
dX §u = 1 ((■ +0uk )= [°u]{pj,
где
5[u ]=[5uj 5wj S3 j Sum Swm S3m 5uk 5wk S3k ].
Аналогично
б 2V
бх2
= [Рм ]{} = 4 [Р1 Р2 Рз Р4 Р5 Рб Р7 Р8 Р9 М
* I
б V 1 N 1( )
= -¡г [8“]{Р«}.
В этих формулах Р1 = Р4 = Рб = Р7 = 0, Р2 = -9б^2 +12^ + 8, Р3 = I(- 24^2 + б^ +1), Р5 = 192Е,2 -16, р8 =-9б^2 -12£ + 8, Р9 = !г (24^2 + б£ -1).
Приведя, таким образом, все матрицы, входящие в формулу потенциальной энергии, к одной размерности [9*9], выполним операции умножения, сложения и интегрирования. При расположении сомножителей в произведении будем исходить из того, что в результате должна получиться матрица размерностью [1*1], так как энергия (и ее вариация)
- скаляр. После соответствующих преобразований получено следующее выражение вариации потенциальной энергии (здесь и в последующем индекс «/» опущен):
ъи =
Е1 0,5 Г, би б
=
1 - 0,5
к------------Ъи + 2 2
бх бх бх2 бх
EJ 0,5
б5 = -гг 1{и[Ъи]{)[Р~]М+[Ъи]{Рм}Рм](и},
I -0,5
к I2 ■
где к = —, /У - радиус инерции сечения относительно нейтральной оси.
В результате сложения двух квадратных матриц и последующего интегрирования получим матрицу жесткости элемента [к. ] . Покажем, например, как найден второй диагональный член матрицы:
Е1 0,5 Е1 0,5 / \7 Е1
к.22 = ~Е3 1 РРб, =-тг I (-9б52 +125 + 8)) = б3,2• -1.
1 - 0,5 1 - 0,5 1
Остальные члены матрицы найдены аналогично. Окончательно матрица жесткости элемента выглядит так:
[к. ] =
Е1
к 0 0 0 0 0 -к 0 0
0 б3,2 18,81 0 -102,4 0 0 39,2 - б^
0 18,81 7,212 0 - 25,б1 0 0 б,8! -1,2!2
- 0,5 0 0 1 0 0 - 0,5 0 0
0 -102,4 - 25,б1 0 204,8 0 0 -102,4 25^
0 1,5 0,251 0 0 I 0 -1,5 0,25l
- к 0 0 0 0 0 к 0 0
0 39,2 б,81 0 -102,4 0 0 б3,2 - 18^
0 - б. 81 -1,212 0 25,б1 0 0 -18,8/ 2 сч 7
(8)
Заметим, что все числа, входящие в матрицу, - точные, без округлений.
Из элементарных матриц формируется глобальная матрица жесткости - таблица коэффициентов левой части системы уравнений. Правую часть системы составляет столбец внешних сил (и моментов), приложенных в узловых точках, к которым нужно
У
3
I
присоединить и распределенную нагрузку. Нагрузку Ц, равномерно распределенную по длине элемента, нужно разделить между узлами элемента. Подставим в интеграл, выражающий вариацию работы распределенной нагрузки, полученное аппроксимационное уравнение прогиба (3):
0,5/ 0,5 , , , ,
I = Ц/ |[ (-8^4 + 2^3 + 4£-1,5)+ )(-2^4 +^3 + 0,5^2 -0,25^)+
- 0,5/ - 0,5
+ 5мт (16^4 - 8Е,2 + 1)+5^к (- 8^4 - 2^3 + 4^ +1,5)+/53] (2^4 +^3 - 0,5^2 - 0,25^)м =
= — ц/ ■5wi + —ц/2 -53, + — • ц/ + — ц/ ■5wk - —ц/2 -53к.
30 ] 60 ] 15 30 к 60 к
Значит, равномерная нагрузка распределяется по узлам элемента таким образом: в крайних узлах прикладываются сосредоточенные силы по 7 ц/ / 30 и сосредоточенные моменты + Ц/2 / 60 , а в среднем узле - только сосредоточенная сила 8ц/ /15 .
Далее приводится полученная тем же методом матрица жесткости элемента пятого порядка. В этом случае появляется множитель 7 в знаменателе и, в отличие от элемента четвертого порядка, некоторые коэффициенты матрицы не изображаются конечными десятичными дробями. Поэтому представим их в виде обыкновенных дробей:
[*. 1 =
Ей
к 0 0 0 0 0 - к 0 0
0 145— 35 18 32—/ 35 0 -102,4 546 / 7 0 СП 1 ^ т | со 03 4 - 32 6 — / 35
0 18 32 — / 35 917 / 2 35 0 - 25,6/ 91 /2 7 0 32 - 6— / 35 1-3./2 35
- 0,5 0 0 1 0 0 - 0,5 0 0
0 -102,4 - 25,6/ 0 204,8 0 0 -102,4 25,6/
0 546 / 7 91 /2 7 0 0 364 /2 7 0 546 / 7 91 /2 7
- к 0 0 0 0 0 к 0 0
0 СО 1 ^ т I со 03 4 - 32 - 6 — / 35 0 -102,4 - 546 / 7 0 14517 35 18 - 32 — 35
0 32 6—/ 35 1-1 / 2 35 0 25,6/ 9 — /2 14 0 18 - 32—/ 35 9 — /2 35
Равномерная нагрузка распределяется по узлам элемента так же, как и в случае элемента четвертого порядка.
При линейном изменении нагрузки по длине элемента
Ц = Ц] (0,5 -£) + Цк (0,5 + £) представим ее как сумму постоянной и переменной частей
Ц = (Ц; + Цк V2 + (Цк - Ц] X,
для каждой по отдельности рассмотрим вариацию работы и найдем следующий закон распределения нагрузки между узлами.
3
/
/
Постоянная часть нагрузки делится между узлами по прежнему правилу, причем для элементов обоих типов одинаково: к j -му узлу прикладываются сосредоточенная сила
7(q, + qk )l/ 60 и сосредоточенный момент (qy + qk )l2 /120, к k -му узлу - такая же сила и
противоположный момент, к среднему узлу m - только сосредоточенная сила 4(q, + qk )l/15.
Переменная часть нагрузки распределяется так:
1) элемент 4-го порядка - к j -му узлу прикладываются сосредоточенные сила -fe - qj )l/10 и момент -( - qj )l2 /120, к k -му узлу -сила (qk - qj )l/10 и момент
-(qk - qj )l2 /120;
2) элемент 5-го порядка - к j-му узлу прикладываются сила -(qk - qj )l /14 и момент
-(qk - qj )l2 / 280, к m -му узлу - момент - 9(qk - qj )l2 /140, к k -му узлу - сила (qk - qj )l /14 и момент - (qk - qj )l2 / 280.
Следует заметить, что в случае распределенной нагрузки переменной интенсивности лучше использовать элемент пятого порядка, который абсолютно точно описывает линейно-переменную распределенную нагрузку.
Предлагаемые элементы можно использовать при расчете плоских рам, состоящих из прямолинейных участков и нагруженных в главной плоскости. В местах сочленений, в тех сечениях, где приложены сосредоточенные силы и моменты, в местах начала и конца распределенной нагрузки образуются узлы, таким образом, каждый участок представляется одним элементом и число элементов получается действительно конечным.
На основе рассмотренных элементов составлена программа расчета плоских рам на алгоритмическом языке «Фортран-90». Программа универсальна, она имеет неизменяемую основную часть, где выполняется решение системы уравнений по определению узловых перемещений (методом Гаусса), определение внутренних силовых факторов и вывод на печать результатов расчета; программа легко перенастраивается благодаря наличию мобильной изменяемой части, в которую входят разные варианты модулей ввода исходных данных (схема, размеры, нагрузки, нумерация узлов и элементов), подпрограммы составления матриц жесткости элементов разного порядка, различные подпрограммы матричной алгебры.
Ниже приводятся результаты расчета рамы (рис. 3) при следующих данных:
l = 2 м, F = 18 кН, q =10 кН/м, M = 16 кНм. На рисунке показана нумерация узлов и элементов (в кружочках со стрелками - указателями направления локальной оси x ). Она, а также геометрические размеры, значения нагрузок, местная нумерация узлов в каждом элементе, углы наклона элементов вводятся с помощью операторов data. Рама состоит из 9 элементов и имеет 17 узлов, значит, число неизвестных перемещений и уравнений равно 51. Углы наклона элементов а отсчитываются от обычного положительного горизонтального направления против хода часовой стрелки, например, угол наклона 1-го элемента
1 -л / 2, 3-го 0 • л/ 2, 8-го 3 -л / 2 и т.д. Каждый элемент имеет двутавровое сечение (№16) с изгибной жесткостью EJ = 1746 кНм и радиусом инерции iy = 0,0657 м.
После ввода исходных данных начинается формирование матрицы жесткости. Перемещения сквозным образом нумеруются в таком порядке:
= 3(1) и4 = и$, " = ,Л/(1) " = 3(1)
и5 =
т >
u3 =&т,
Сначала матрица 51*51 обнуляется, затем, в порядке нумерации, вызываются элементарные матрицы жесткости. Вычисленная матрица жесткости элемента приводится к виду, соответствующему ориентации элемента в глобальной системе:
[*. ]=[с ]-1[*. ][С ],
cos a sin а 0
где
- sin а cos а 0
матрица перехода,
0 0 1 1
[к. 1- приведенная (к глобальному виду) элементарная матрица. Здесь показана только одна из трех одинаковых диагональных подматриц, остальные элементы матрицы [С ] размерами 9 х 9 равны нулю.
После вызова и описанной обработки первая матрица занимает свое место в глобальной матрице - левый верхний угол из 9 х 9 = 81 ячейки. Члены второго элемента расположатся вслед за первым, вдоль главной диагонали глобальной матрицы, «зацепившись» за последний тройной квадрат первого элемента своим первым и т.д. В третьем узле сходится три стержня, и в ячейки общего тройного квадрата попадут коэффициенты от трех элементов. Содержимое общих для нескольких элементов ячеек складывается.
После подготовки матрицы жесткости формируется вектор правой части системы уравнений, в соответствии с местом приложения и направлением внешние силы занимают свои места в столбце правой части системы. Вводятся граничные условия. В данной задаче жестко защемлены узлы 1 и 17, поэтому
В программе каждое граничное условие реализуется так. Возьмем, например, условие и50 = 0. Все коэффициенты 50-й строки и 50-го столбца матрицы и 50-й член столбца
правой части полагаются равными нулю, а коэффициент в ячейке на их пересечении -единице.
Подготовленная таким образом система уравнений решается методом Гаусса, и находятся перемещения узлов. Для определения внутренних силовых факторов нужно перейти от глобальных перемещений к локальным. Это делается с помощью упомянутой матрицы перехода:
М=[с ]{и},
где {и} - вектор локальных перемещений,
(и } - вектор глобальных перемещений.
Внутренние силы находятся по следующей матричной формуле:
N к 0 0 0 0 0 к 0 0 "
Оі 0 108 30/ 0 -192 0 0 84 -18/
02 0 8 - -18/ 0 192 0 0 -108 30/
Мо » 0 сч сч - -8/2 0 32/ 0 0 -10/ 2/2
Мс 0 81 /2 0 -16/ 0 0 6/ - /2
М к. 0 -10/ - 2/2 0 32/ 0 0 - 22/ 8/2
ип
и„
где N - нормальная сила (постоянная по длине элемента),
01 - поперечная сила в начале элемента,
02 - поперечная сила в конце элемента,
М0 - изгибающий момент в начале элемента,
Мс - изгибающий момент в середине элемента,
Мк - изгибающий момент в конце элемента.
Результаты расчета представлены на рис.4 эпюрой изгибающего момента, построенной на сжатых волокнах.
Рис. 4 - Эпюра изгибающих моментов, кНм Литература
1. Чирас, А.А. Строительная механика/ А.А.Чирас. - М.: Стройиздат, 1989. - 256 с.
2. Сегерлинд, Л. Применение метода конечных элементов/ Л.Сегерлинд. - М.: Мир, 1979. - 392 с
© А. Х. Валиуллин - канд. техн. наук, доц., проф. каф. теоретической механики и сопротивления материалов КГТУ, 1ш8ш&к81и.т.
