Криволинейный балочный конечный элемент, учитывающий геометрически нелинейные деформации Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАпИСКИ ЦМИ

Том XXII

1991

N5

УДК 629.7.015.4.023

КРИ^ШИНЕЙИЫЙ БАЛОЧНЫЙ КОНЕЧНЫЙ ЭЛЕМЕНТ, УЧИТЫВАЮЩИЙ ГЕОМЕТРИЧЕСКИ НЕЛИНЕЙНЫЕ ДЕФОРМАЦИИ

С. Н. Зайцев

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

В современных конструкциях широкое применение находят балочные конструктивные элементы. В ряде случаев при проектировании таких элементов необходимо учитывать геометрически нелинейные деформации. К таким конструктивным элементам относятся, например, лопасти винтовентиляторных дваигателей и вертолетов, подкрепляющие ребра оболочек. Моделирование этих элементов наиболее естественно осуществлять с помощью балочных конечных элементов. Разработке таких конечных элементов посвящено большое число работ. Однако многие предлагаемые в публикациях балочные конечные элементы не обладают достаточной универсальностью и не пригодны 1} случае произвольной пространственной геометрии, произвольной степени анизотропии, формы поперечных сечений, конечных прогибах и углах поворота.

Часто применяется упрощение, связанное с линеаризацией перемещений элемента по отношению к приращениям узловых углов поворота, т. е. углы поворота считаются бесконечно малыми на каждом шаге увеличения внешней нагрузки [3]. Такое упрощение имеет, однако, тот недостаток, что приводит к более грубой геометрической матрице жесткости элемента и, как следствие, существенному увеличению общего числа итераций при нелинейном расчете.

Авторы работ [1] и [2] предложили сходные формулировки криволинейного геометрически нелинейного балочного конечного элемента, в которых строго учитываются конечные прогибы и углы поворота, принимаются во внимание деформации поперечного сдвига и которые могут применяться для моделирования пространственных балок, сделанных из

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

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

1. Балочная модель геометрически нелинейного деформирования сплошных тел. Условимся отсчитывать конечные перемещения точек тела от некоторого известного (не обязательно равновесного) деформированного состояния тела, помеченного индексом t. Это условие связано с выбором модифицированного подхода Лагранжа [3^ 4] при рассмотрении процесса геометрически нелинейного деформирования. Применение модифицированного подхода позволяет получить в данном случае более простую конечно-элементную формулировку. Однако приведенная в работе математическая и конечно-элементная формулировка может быть легко модифицирована для применения общего подхода Лагранжа [3, 4]. При общем подходе отличие состоит только в том, что роль состояния / играет начальное недеформированное состояние, не изменяющееся в процессе итераций.

Считаем, что в состоянии / фиксирована криволинейная система координат 61, 62, 6з (рис. 1).

Введем обозначения:

Я = Я' (61, £2, 4з) — радиус-вектор произвольной точки балки в состоянии

? = $'(61 О О) — радиус-вектор произвольной точки упругой линии балки в состоянии / (выбор упругой линии определяется соображениями наилучшего соответствия принимаемых балочных гипотез напряженного состояния условиям реального деформирования);

ёг(60, ёз(60 — единичные, взаимноортогональные векторы криволинейного базиса вдоль упругой линии балки и лежащие в сечениях балки, не являющихся в общем случае в точности перпендикулярными к упругой линии из-за наличия малых деформаций поперечного сдвига в состоянии / (в первоначальном недеформированном состоянии сечения задаются ортогональными) :

(1)

Введем кинематическую гипотезу плоских сечений, учитывающую влияние сдвига [5], согласно которой радиус-вектор Я и векторы криволинейного базиса любой точки тела в начальном состоянии / и

,1 >

1 = (ёгёу.

Рис. 1

произвольном последующем деформированном состоянии, помеченном индексом t + ДЛ имеют вид:

2 ' с2 "г &3 • ег

? + е2 + • /Г" « + Е,. + Ь.

ё\ = = ех + £2 • е2,| + £з • ¿ал I

Я* п . о Я'+М -»+Л'

8: =(а = е*, а =2 3, 8а = = ; (<?2 .^' = (е3- е-)' =1, (<?2*«Га), = О;

(<?2-е2) =(е3.е3) =1, (е2 • ез) =О.

\> (2)

Считается, что при деформациях векторы а = 2 3 остаются орто-нормированными и поворачиваются как целое на конечный угол, отсчитываемый от состояния Поле перемещений балки имеет вид:

^Д* == Г<+Д< — г1, ы'+4< = — Й'= +

+ ь + - щ (3)

и, следовательно,

Ц+м = Г'+м = ё + (4)

где — перемещения ■ точек упругой линии балки.

^обЩенными перемещениями балочной модели могут считаться, таким образом, одномерное поле перемещений и одномерные векторные поля (ё+4' — ё!?) и — Однако векторы ё!?, ^ не могут перемещаться в пространстве при деформациях независимо друг от друга. Совместность деформирования векторов ё обеспечивается введением углов конечного поворота в точках упругой линии. Кинематика этих углов будет рассмотрена ниже при описании конечно-элементной аппроксимации геометрии и поля перемещений элемента. Для последующего важно пока

к(5)

лишь то, что векторы £4+4' являются нелинейными вектор-функциями трех аргументов, имеющих смысл компонент углов конечного поворота.

Тензор конечных деформаиий Грина сплошного тела [3, 4]

в состоянии t + Ы для заданной кинематической гипотезы (2) и в случае «тонкой» балки, дЛя которой отношение характерного размера сечения А к радиусу кри&изны # п/Н ^ 1 (ранное предположение позволяет пренебречь в выражении для е^4' членами, содержащими У н £2, £з), имеет вид:

/+д/ 1 (-*+л< -А

е'+А< - _|г. + бз • К2+М, 2е{+*' _ у'2+&1 - Ъ • г'+л'. 2е{з"л< = + и ■ г'2Г' - ' - е'+л' = О ,

=±(?г • ^ - г;. ,1). =- .?<,

« • ¿Г -1 ■ е3, Г'+л' = е3+&' • - ¿<3...

1Г1+АI -<+А' . -<

Обозначим

е/н/ = Ие'+М11 = 1е'\'2'\'зТК2Кз1;+м,

[6Х 1J — вектор обобщенных деформаций балочной модели. Везде далее будем предполагать деформации малыми и пренебрегать изменениями объема, площади и длины при деформациях.

Для определения обобщенных напряжений балки рассмотрим виртуальную работу внутренних напряжений в состоянии t + Ы. Используем 2-й тензор напряжений Пиолы — Кирхгофа и тензор деформаций Грина

[3, 4J:

би'+"= ( а^-бе^.^, I, /= 1 +3. (6)

Здесь и в дальнейшем по одинаковым верхним и нижним индексам подразумевается суммирование. Используя (5), имеем:

би+м= ^ (¿1. • Со = Н^'бе + Р26'\'2 + ^Збуз +

0(+Д/ 0

М'бТ + М2бК2 + МЗбКз)'+Л'с/ = ((а1 • бе^'+Л<. а,

о

^(+4/ = с18 • М, dS = сС12 • сС£з ,

с=/сб„

Р'+Д( = д/ * , Р/+Д/ =,$СТ'+Д< • dS, Р?+д< = 5а,+д/ • dS,

5

>(7)

М;+м = ^ (о1 % — ° %) + • М2+М — \

= — Ц+М^

Выполненное в (7) преобразование интеграла по объему в интеграл по длине балки справедливо, вообще говоря, лишь для прямолинейных балок постоянного сечения. Однако если предполагать, что балка «тонкая» (в указанном выше смысле) и сечение изменяется достаточно постепенно вдоль балки, то можно пренебречь погрешностью этоги преобразования. Данное замечание относится также к последующим результатам, при вы-водекоторых используется преобразование объемного интеграла в интеграл по длине. Обозначим

Ч + Ы

=

/ + А<!

[6Х 1] — вектор обобщенных напряжений балочной модели.

При переходе к новой координате ^ вдоль упругой линии векторы ^общенных деформаций и напряжений преобразуются, как это следует из (5), (7) и тензорного характера величин о7+д, и е/+Л', следующим образом:

Е;+м — В * 8/+м. °(+д/ = В 1 •

В

О

ь>

Считаем, что векторы обобщенных напряжений и обобщенных деформаций связаны линейной зависимостью — обобщенным законом Гука:

(+Д/

а', + Нц. е/+Л', (,/ = 1 Ч- 6.

(9)

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

Упругие постоянные матрицы Н находятся отдельно для каждого сечения тела и полагаются равными упругим постоянным эквивалентной прямолинейной, однородной по длине балки, имеющей такое же сечение. Жесткостные свойства этой эквивалентной балки определяются с учетом нестесненной депланации сечений, а также принятия обычных балочных гипотез о характере напряженного состояния и зависят только от конструктивных и геометрических параметров сечения. балки [6].

Таким образом, при построении балочной модели считается, что в каждом сечении криволинейной балки переменного сечения накопление упругой энергии при деформировании происходит так же, как и для эквивалентной прямолинейной, однородно^балки.

В соответствии с данной гипотезой вычисление деформаций тела на основе (5) должно быть дополнено, чтобы учесть влияние депланации [6].

Компоненты матрицы Н зависят от выбора координаты £1 вдоль упругой линии балки. Закон преобразования матрицы Н при переходе к новой координате можно получить,используя инвариантность произведения отбе и (8):

о'тбе' = е'тН'бе' = етВтН'Вбе,

а'т6в' — отбе = етНбе, Н = ВТН'В.

О

ь

ь

ь

Матрица Н', определяемая при рассмотрении деформаций прямолинейной, однородной по длине балки, соответствует случаю, когда роль координаты 61 играет длина дуги упругой линии. Поэтому матрицу Н, входящую в (9), можно получить с помощью (1О) и (8), учитываял что:

' (11)

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

ба'+д' = б ^ + 62 бё+Ы + 6збЙ+ы, (12)

6Й+4' = бфХё+ы, а =2 -т 3. (13)

Согласно (5), имеем:

бе'+4' - ¿Г' • б Г „ - . + еТ4' ■

1

= ^^ + ¿Г • ,, 6Т'+А' = еТ'■ 6*1Г + С• «¿Г. (И)

= + еЗГв*?,,. _ ¿Г'бС —Г б #,.]

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

5 ((1681) 'аV - 5 (ТбД'+М аи - 5 (ГбД'+М аэ = ¡¡' [ 1б +

"<+4( "'+41 5[+А/ 5Г

1

+ Ябч», + (РХ ё.) 6сР] 'Н' ¿1,- 5 (РбФ+ЕбсР)'Н' ¿6, _

-(Рб ^ + Сбф)

!Г

= О ,

-(Рб* -Сб^)'Н'

5Г

Р= • + /*в2 + + М2езл - М3ё2..). (15)

+ ^а)'), + 6;) х

X fdS + + 6;) X ^Оу) ,

р±= 5 тм, = 5 (ьл+6;) х Tds. J

Осуществляя процедуру варьирования, получим обобщенные нелинейные уравнення равновесия в состоянии t + Д^

1,1+ [5=0, (16) = ± Р±, = ± .

Полученные уравнения для ? и М совпадают с известными нелинейными уравнениями равновесия криволинейной балки [5], которым удовлетворяют физические силы и моменты в сечениях деформированного тела. Поэтому величины Г и М совпадают с физическими силами и моментами в сечениях. Следствием этого являются следующие равенства, выражающие связь.между обобщенными напряжениями Р, М' и внутренними усилиями Р, М' в сечениях:

?= Я ал/ей" + Р2-ё2 + Р3-ё3, # = лг'-а/УёТГ + м2-ё2 + м3-ё3,

Р = е,,(Г1 + М2(ёз,.ё.)/е.. - М3(ё2,\ё\ ) /ем), ( ,

р2 = (Г2 + м2(ёз,.ё2)), Р3 = -{ё.-) 1 М1 , м2 = еп-м2, М* = е„.М3 .

Установленное различие обобщенных напряжений балочной модели и физических" сил и моментов в сечениях является следствием криволи-нейности балки и известного факта, что применение метода перемещений влечет за собой выполнение уравнений равновесия только в некотором интегральном смысле, в «среднем:., в то время как поле напряжений в теле оказывается локально неравновесным.

Сформулируем теперь в обобщенных переменных принцип виртуальных работ для балочной модели нелинейного деформирования тела. Этот принцип является основой для последующего построения шагово-итерационного процесса расчета нелинейного деформирования тела. Принцип записывается для равновесного состояния тела t + Д/, но, в отличие от приведенной ранее формулировки этого же принципа, переме-щеия и углы конечного поворота отсчитываются от предыдущего известного деформированного состояния ^ Используя (7) и (12), получим (для краткости не включаем усилия на торцах):

( (а''.6е17)'+Л'^- 5 (Г-бйГ"^- ^ (Т■ &й)'+Ы(¡Б =

°(+4< "/+4/

= 5(0' - 6е,)'+Л^/ - $(Р. 6 # + С2 • б- + ёз- = 0, (18)

о о

где

р+м= + ё'Г= 6 + а= 2 З. (19)

Заметим, что вариации 6ё,+А'уже нельзя записать в виде 6фХ£&+Л', поскольку векторы нелинейно зависят от углов конечного поворота сечений балки, которые отсчитываютсй от состояния По этой причине среди обобщенных внешних усилий вместо привычных моментов внешних сил появляются силовые факторы С"'.

2. Инкрементальная теория расчета геометрически иелииеАиых де-

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

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

Линеаризованные уравнения равновесия шагово-итерационного процесса получаются с помощью инкрементального вариационного принципа, построенного на использовании общего или модифицированного подхода Лагранжа [3, 4]. Сформулируем этот принцип для балочной модели деформирования. Будем отсчитывать перемещения и деформации балки от известного, деформированного состояния, помечаемого индексом /. При использовании общего подхода Лагранжа этим состоянием является всегда начальное недеформированное состояние. В случае модифицированного подхода . Лагранжа в качестве состояния / может быть выбрано любое известное деформированное состояние, достаточно близкое к искомому равновесному состоянию на данном шаге внешней нагрузки, например равновесное состояние, соответствующее предыдущему уровню внешней нагрузки. При выводе инкрементального принципа в качестве переменных используются 2-й тензор напряжений Кирхгофа и тензор деформаций Грина.

Искомое равновесное состояние для заданного нового приращения нагрузки помечаем индексом / + об/. Для получения инкрементального вариационного принципа запишем принцип виртуальной работы балки (18) для состояния / + М (для краткости не включаем усилия на торцах):

¡(о1 • 6е,)'+А'<И- \{Р. Ш + + бё3)'+4'<*/ = 0.

о о

После чего, все переменные, зависящие от компонент обобщенных перемещений, разложим в ряд в окрестности известного состояния с индексом /л по степеням приращений компонент обобщенных перемещений от состояния /„ до / + об/ (/„ — неравновесное состояние, найденное в реэулотате предыдущей п-й итерации при итерациях на данном шаге приращения внешней нагрузки; /о совпадает с /^л л

В ' приводнмых ниже выражениях символы 6, 62, ... означают соответственно линейные по приращениям компонент обобщенных перемещений, квадратичные и т. д. члены разложений. учитывается также нелинейная зависимость векторов ё« от углов конечного поворота. Считается, что внешние нагрузки могут зависеть от перемещений:

б8,'+м=б( в;-+ее;-+82е'г+...) = б( ее;-)+б( л;-)+... , б(ве;-) = бе;» ;■

<К+А< = 4 + НЦе'+м = а' + Нце'.; + + ... = о,. + И • ... , '

(о . беУ+М = О}., бе,'" + б[ +(Я"^" + 20;.Л,'")] + ... ,

=б( + = 6Д^ = Л(20)

б-г=б( +«¿чбй+...) = ++..., = +&'+..., сг= Й++...,

+С2бё2+Сзбё^ 'Н'=+С2бё2+Сзбёз)'" +

++§с2б£2+§Сз- бёз)"+б(е/ё2) + б(ёзБёз) + ...

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

б|"$+(Н ■ • + 2а'б2е)'"1 d/ _ ([ + §С2бё2 + аВД,)'" +

'-О О

'пТ 1 I

1;-б(с2б?ё2+СзБёз) ] d/ = ¡(р.+ с2^ бё2 + Сз- бёз) (21)

о

1

-Но' . бе;) " .

о

С помощью варьирования (21) получаются линеаризованные уравнения равновесия относительно приращений обобщенных перемещений, соответствуюшие (п + 1)-й итерации, которые позволяют найти (п + 1)-е приближение по известному п-му.

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

С помощью (5) можно получить выражения для бе,- " и б2Е(", фигурирующие в (21) и необходимые для дальнейшего рассмотрения:

8в ■ = е" • Д *, = (¿1 + Д , = + + (22) $К'2п = (в,8в8-1 + А*, • ¿з..)'". (- в,8?,,,' - А*;, •

§У» = -|(Д= • + • 5^'", (23)

^з" = (Д^.1 • 8?3 + е, • 8%) \ бУ" = (8ё3 • | + ¿3^,, +

'п '«

3. Конечно-элементная формулировка решения нелинейной задачи.

Инкрементальный вариационный принцип (21) позволяет применить стандартную процедуру МКЭ аппроксимации геометрии и поля перемещений и получить линеаризованные матричные уравнения равновесия балки.

Введем в каждом элементе балки следующую аппроксимацию геометрии и поля перемещений для произвольного деформированного состояния [3]:

р+м = ё+м = ё+м, Ф'+м = W/H', (24)

i = 1 n, а. = 2 3,

где N — функции формы, n — число узлов элемента.

Будем также считать, что узловые векторы ё+м получаются в результате конечного поворота соответствующих векторов ^ исходного деформированного состояния t. Угол конечного поворота в узле определим аналогично работе [1] следующим образом. Введем в узле дополнительный ортонормированный базис äj и выполним три следующих друг за другом конечных поворота вокруг осей ä1, а2, аз на угол ч>3

соответственно. Новое положеиие вектора обусловленное этим конечным поворотом, можно получить, если найти новое положение ортов äj, вызванное конечным поворотом, и учесть, что вектор Ц+А' занимает относительно повернутых ортов а}+А1 такое же положение, какое вектор £ занимал относительно старого базиса äj:

£ = E'aa', E'a = \(ё.. а.) (г. а)(г. аз)!', }

а' = \ Щ ЩЩ |', а'+м = ф'+м а| (25)

ё+м = Ё^а'+м = Е'аФ'+м'а1, J

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

ф=

cos92cos93

Sin 9,Sin q^COS ч>з-

— СОЭф^Шфз C0SфlSinф2С0Sфз + + SiПфlSiПфз

COSф2SiПфз Sin ч>1 Sin ч>з +

+ coSфlCOSфз COS ч>1 Sin Ч>2 Sin ч>з —

— Sin9.cos^3

— Sin^2 Sinq>iCOS92 COS ч>1 COS Ч>2

• (26)

Дополнительные узловые базисы 5 можно выбирать произвольно. В частности, можно использовать постоянный базис пространственной системы координат (глобальный базис) или изменяющийся базис, орт 51 которого направлен по а остальные два орта лежат в перпендикулярной плоскости и связаны каким-либо способом с ортами ^ (локальный базис).

Аппроксимация (24) — (26) обладает, как можно убедиться, свойством конечного движения как твердого целого. Это качество должно обеспечиваться, так же как и его аналог, в линейном расчете. Невыполнение этого свойства может отрицательно сказаться как на скорости сходимости итерационного процесса, так и на конечном результате.

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

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

А _ Следствием (24)—(26) являются выражения для вариаций б2е::, ДФ, которые используются при формировании основных матриц и векторов балочного конечного элемента на основе вариационного принципа (21). Заметим, что принятые при построении итерационного процесса

А Л- .

соглашения приводят к наиболее простой записи для бё::, б2ёа, которые при другой методике получаются довольно громоздкими:

= N ■ бй = N - в'ы. &ф|". а = N ■ Ё1 Дф{, 1

= N . = N-В<ы. Гф\- . а.: = N £21- . Дф? Дф/ , (2.7) Д Ф = NД Ф, = N-Дш! . а,,«' = 1 п, j, k, / = 1 3, а. = 2 3 .J

Выражения для векторов £2^/,- легко получить при рассмотрении

первой и второй вариаций матрицы конечного поворота Ф в состоянии с учетом того, что углы ^ в состоянии ¿„ являются нулевыми.

Используя введенную МКЭ-аппроксимацию и осуществляя процедуру варьирования согласно выражению (21), получим линеаризованные матричные уравнения равновесия балки, которые имеют следующий вид:

(А^+Ка+А^-Ди(п+1) = /5(п)-5("), (28)

где К(£П) — линейная матрица жесткости в состоянии ¿я; Ю„п) — матрица геометрической жесткости в состоянии ¿я; Л^Г' — матрица жесткости внешних сил, зависящих от перемещений. В румме эти три матрицы составляют матрицу касательнои жесткости; Р — вектор внешних узловых сил в состоянии ¿„; — вектор внутренних узловых сил в состоянии

Инкрементальный вариационный принцип (21) является также основой для получения линейных уравнений малых колебаний в окрестности положения статического равновесия и с учетом дополнительной геометрической жесткости от действия напряжений равновесного состояния. В этом случае во внешнюю нагрузку должна войти дополнитэдьная инерционная составляющая. Кроме того, состояние tn везде в (21) должно быть заменено на состояние t, поскольку линеаризация осуществляется относительно состояния статического равновесия. Получаемые с помощью (21) матричные уравнения колебаний МКЭ записываются аналогично (28):

ММ + СМ + (Кь Ч Ка + КР)М = АР, (29)

где М — матрица масс; С — матрица демпфирования от действия доба' вочных внешних сил, зависящих от скоростей движения точек тела; ДР — вектор внешних узловых сил от действия дополнительной возмущающей нагрузки, не зависящей от перемещений.

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

Линейная матрица жесткости элемента формируется на основе (22), (24), (27), (10): .

ие= |Дю1Дг0?Дш3Дф!Дф^ (30)

где ые — [6п Х 1] вектор узловых перемещений элемента, п — число узлов элемента.

При формировании геометрической матрицы жесткости используются (2З), (24), (27), (9):

(31)

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

С помощью тех же соотношений получаются другие матрицы и векторы элемента:

- ^ + [;ё2бё2 + [;ёз&з) + б( ё2б2ё2+ёзбё)"] ш=бит ке„ие\

о

^ '

5 + С2б£2 + ёзбёз)" Ш = 6«; . (32)

о

5 ( О' •бе,) "Ш = б и] •Бе.

Построение матрицы масс элемента рассмотрим более подробно:

"5, Са =

S 5

= 5-рДы "5, С. =5-р^Ды <1Б,

дГ+Л<=Ш'+гМ+1М+...,р'=б'2=с'3=0,

$р< = _ Г^(рЛ8 + . 5 р + ^«з • 5 Р • 6з<«1 ,

I .5 5 Я -л

ЕС = - [Др+ |р, |2- + ^р- !з- б.^], - 5е [ (£р- б г + . £С2 . бё2+^ёз • бёз) '+б( ё282ё2 + ёз$2ё,)'ш] _

о

= 6итеМеДйе.

Ш)

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

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

В случае использования общего подхода Лагранжа описанная процедура построения элемента и итерационного процесса получает незначительные изменения, связанные главным образом с более сложными

выражениями для векторов £1!";, Ё2Мми в выражениях (27), поскольку усложняются первая и вторая вариации матрицы конечного поворота Ф (26) в состоянии /„, так как углы ^ в этом состоянии становятся отличными от нуля.

4. Изгиб консольной балки. В качестве первого примера использования предложенного балочного конечного элемента рассмотрим задачу об изгибе консольной балки под действием изгибающего момента на конце, как показано на рис. 2. Конечно-элементная модель балки состоит из четырех трехузловых элементов. Данный пример взят из работы [1]. Согласно теоретическому решению балка деформируется по дуге окружности, причем величине момента М = {(л£//1) соответствует угол сектора.равный = {л. Полный момент, прикладываемый за ряд шагов, соответствует коэффициенту { = 1,8 и углу 1,8л.

Линеаризованное уравнение равновесия данной задачи, решаемое на каждой итерации, имеет вид:

где вектор внешней нагрузки Р имеет единственную отличную от нуля компоненту, равную величине изгибающего момента на данном шаге и соответствующую углу поворота вокруг оси У в последнем узле балки.

Так же как ив [1], равновесие считалось достигнутым, когда максимальная компонента вектора невязки в правой части (34), отнесенная к текущей величине момента, становилась меньше 0,001.

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

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

Таблица 1

Парамер f б, 6, (теория) бх бх (теория) Невязка А Число итераций

0,15 -2,776 -2,776 -0,441 -0,439 9,3Е—6 6

0,3 -5,253 -5,249 -1,710 -1,700 1,0Е—5 6

0,45 -7,169 -7,160 -3,649 -3,616 1,5Е—5 6

0,6 -8,336 -8,333 -6,020 — 5,946 • 2,6Е—5 6

0,75 -8,665 -8,694 -8,532 — 8,400 4,9Е—5 6

0,9 -8,174 -8,281 -10,88 —10,"69 9,0Е—5 6

1,05 -6,995 -7,231 -12,80 —12,57 1,6Е—4 6

1,2 -5,349 — 5,758. -14,07 —13,87 2,6Е—4 6

1,35 -3,528 -4,114 -14,58 -14,52 3,7Е—4 6

1,5 -1,849 -2,546 -14,34 -14,55 4,7Е—4 6

1,65 -0,607 -1,264 -13,50 -14,06 7,3Е—4 6

1,8 — 0,027 -0,405 -12,33 -13,25 3,7Е—4 7

Иг >го: 73

Таблица 2

Результаты работы | 1 |

Параметр f б, б, Невязка А Число итераций

О,2 -3,623 -0,765 7,9Е—5 7

0,4 -6,574 -2,899 1,2Е—4 7

0,6 -8,328 • -5,955 2,6Е—4 7

0,8 -8,620 -9,293 5.5Е—4 7

1,0 —7,507 -12,22 3.5Е—3 7

1,2 -5,376 — 14,12 30Е-9 8

1,4 -2,877 -14,66 3,2Е—7 8

1,6 -0,797 — 13,85 7,9Е—5 8

1,8 +0,123 -12,18 1,0Е—6 9

Итого: 68

5. Потеря устойчивости криволинейной полосы. Рассмотрим тонкую криволинейную балку, упругая линия которой совпадает с дугой окружности радиуса Я и углом раствора р (рис. 3).

Известно теоретическое решение для критической величины изгибающих моментов, приложенных на торцах балки, под действием которых происходит боковая потеря устойчивости балки [5]:

М,,2 = [(£/ + GJ) ± -^Ё/—.

Эта задача была решена с помощью МКЭ, и проведено сравнение с теоретическим решением. Для моделирования балки использовался балочный конечный элемент, построенный в работе.

Воспользуемся бифуркационным критерием потери устойчивости, согласно которому при критической величине действующей нагрузки наряду с основным невозмущенным состоянием равновесия t существует смежное бесконечно близкое возмущенное равновесное состояние t + А/. Применим для этих двух состояний инкрементальный вариационный принцип (21). Учтем, что роль состояния tn в данном случае будет играть равновесное состояние .а также то, что правая часть в записи этого вариационного принципа представляет собой сумму виртуальных работ внешних и внутренних сил в равновесном состоянии t и, следовательно, в соответствии с принципом виртуальных работ обращается в ноль. Вводя балочную МКЭ-аппроксимацию, в которой узловые перемещения относятся к локальным базисам получим соответствующую вариационному

принципу линеаризованную матричную задачу на собственные значения для определения критических сил и форм потери устойчивости:

(К + ЛК> = О.

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

рх = ? = р3 = М' = М2 = 0, м3 = м, Я = ^з=м' = м2=0, М3 = А = - = -Ц-.

еи еи Я

Матрица Ка балочного элемента формируется для известных М\ М2, Мз так, как это было описано ранее.

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

В табл. 3 приведены критические величины моментов М 1,2, полученные в случае разбиения балки на два четырехузловых элемента, и сравнение с теоретическими значениями для различных углов 6.

Таблица 3

Угол 11 М. М| (теория) М2 М2 (теория)

45 90 135 180 +373,41 +225,94 + 177,10 + 153,13 +374,28 +225,87 + 176,69 + 152,31 —221,77 —"73,58 —24,42 +0,003 —221,97 -73,56 -24,38 ■ О

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

ЛИТЕРАТУРА

1. S u r а n а К. S., S о r е m R. М. Geometrikally поп-Ппеаг formulation for three dimensional curved beam elements with large rotations.— Int. J. Numer. Meth. Eng., 1989, vol.28.

2. D v о r k.i n Е. N., О n а t е Е., OI i v е r J. On non-linear formulation

rotation increments.- Int. J. Numer. Meth. Eng., 1988, vol. 26, N 7, 1597-1614.

3. В a t h e К. J. Finite element procedures in engineering analysis.— Prentice-Hall, 1982.

4. В а с и д з у К. Вариационные методы и' теории упругости и пластичности.— М.: Мир, 1987.

5. W е m р n е г G. Mechanics of solids. McGraw-Hill, 1973.

6. П о л и щ у к В., Ч у б а н ь В. Уравнения упругой деформации балки с учетом нестесненной депланации в форме метода конечных элементов.— Ученые записки ЦАГИ, 1984, т. 15, N2 1.

Рукопись поступила /7//V /990 г.