Применение метода конечных элементов к расчету больших перемещений плоской линейно-упругой конструкции Текст научной статьи по специальности «Физика»

УДК 531

ПРИМЕНЕНИЕ МЕТОДА КОНЕЧН Ь I X ЭЛЕМЕНТОВ К РАСЧЕТУ БОЛЬШИХ ПЕРЕМЕЩЕНИЙ ПЛОСКОЙ ЛИНЕЙНО-УПРУГОЙ КОНСТРУКЦИИ

А.О. Щербакова1

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

Ключевые слова: геометрическая нелинейность, большие перемещения, конечные деформации, чистый изгиб, плоское напряженно-деформированное состояние.

Введение

Современная механика сплошной среды отличается большим разнообразием подходов к решению задачи о расчете больших перемещений конструкций, разнообразен также и используемый тензорный арсенал. Для описания напряженного состояния материала наиболее широко используется тензор напряжений Коши [1—3], наследуемый практически без изменения из геометрически линейного подхода, а также тензоры Пиолы-Кирхгофа первого и второго рода [2, 3]. Еще большее количество тензоров описывает деформированное состояние материала, например [1-5]. Сюда относятся тензоры конечных деформаций Коши-Грина, Грина-Лагранжа, Альманси и др. Кроме того в литературе, а также в некоторых пакетах вычислительных программ (в частности, в пакете ЛК8У8), для расчета конечных деформаций используют тензор логарифмической деформации (тензор Генки) [5, 6], что делается, как показано в работе [7], не всегда корректно. Подобное обилие мер напряжений и деформаций существенно усложняет решение, причем даже для линейно-упругих конструкций. Поэтому целью данной работы является создание такой модели расчета конечных деформаций плоской конструкции, которая основывается на тензорах, наиболее понятных с точки зрения физического смысла их координат.

В данной работе статическая сторона задачи решается на основе тензора напряжений Коши. Модель деформационных свойств линейно-упругого материала наследуется из геометрически линейного подхода и описывается законом Гука. Для описания геометрической стороны задачи использован тензор деформации, представляющий разность тензора растяжения и единичного тензора, где тензор растяжения определяется с помощью полярного разложения [8, 9] тензора дисторсии.

Описание расчетной модели

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

1 Щербакова Алла Олеговна - кандидат технических наук, кафедра прикладной механики, динамики и прочности машин, Южно-Уральскийгосуда£ственный_униве£ситет.£-тш1:^11а8сЬегЪакоуа@Цй.т^^^^^^^^^^^^^^^^^_^^^^^^^^^^^^^^^^^_

Геометрические соотношения модели. Для рассматриваемой геометрически нелинейной задачи одним из ключевых моментов является задание линейной связи между дисторсией и смещениями узлов конструкции. Однородность деформированного состояния элемента означает наличие линейной связи между начальными векторами волокон 10 этого элемента и векторами его волокон I в деформированном состоянии:

I = ¥-1»- (1)

Оператором этой связи является двухвалентный тензор дисторсии Р. В плоской задаче он представляет сумму четырех диад

¥ = ^1^1^1 + ^2^1^2 + ^з^2^1 + ре 2е 2, где числа Рь Р2, Рз, Р4 - координаты тензора дисторсии в декартовом базисе {е1, е2}.

Тензор дисторсии конечного элемента определяется векторами смещений его вершин (рис. 1). Здесь штриховыми линиями показан элемент ЛВС (заданный векторами а и Ь) в начальном состоянии, сплошными - в деформированном состоянии. Смещения его вершин Л, В и С определяют векторы qЛ, ^ и ^ соответственно. Тензор дисторсии этого элемента определяется с помощью выражения (1), записанного для векторов а и Ь:

¥ • а = а + q В - q Л , ¥ • Ь = Ь + q С - q Л .

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

Рис. 1. Дисторсия конечного элемента

' р' д!4 д Л 1" 2 а 1 2 0 - Ь2 0 а2 0 "

р =[Ь ] дВ + 0 [ь ]= 1 а^ - Ь^ 0 Ь1 0 - а1 0

Рз д 2В 0 , [Ь1] = , а2Ь1 2 Ь1 а1 1 0 Ь2 - а2 0 - Ь2 0 а 2

[ р ] дд 1 1 0 а1 - Ь1 0 Ь1 0 - а1

. (2)

В случае если какое-либо перемещение узла конечного элемента запрещено, мы его исключаем из столбца [д] в выражении (2), удаляя при этом соответствующий столбец матрицы [Ь1]. Например, если для элемента на рис.1 запретить горизонтальное смещение узла В (д1В = 0) и вертикальное смещение узла С (д2С = 0), то его степень свободы, вместо шести, станет равна четы' Ч2Л Ч2В д1С]Т, а матрица [Ь^

0

ц - Ьх 0

рем. Столбец перемещений при этом примет вид [д\Л д2Л д2В д1С]Т, а матрица [Ь^ будет равна

[А ] =

1

а 2 Ь1 - а1Ь2

^2 - а 2

0

0

а1 - Ь1 Ь1

Для конструкции, состоящей из нескольких элементов (их число обозначим к), тензор дисторсии удобно записывать в матричной форме в виде столбца [Р], в котором последовательно в виде блоков перечислены координаты тензоров дисторсии каждого элемента:

Р]=Цр1 Р2 Рз р] Р1 р Рз РаI ••• Р Р Рз РаIГ •

Для конструкции выражение (2) принимает вид

[Р] = [Ь] [д] + [/х]. (3)

Здесь матрица [д] представляет столбец чисел, содержащий смещения узлов конструкции [д] = [д1 д2 ... дт]Т, где число т представляет степень свободы конструкции. Матрица [/1] - это столбец, содержащий блоки координат единичного тензора:

[/1 ] = [ [1 0 0 1] [1 0 0 1] ... [1 0 0 1] ]Т.

Количество этих блоков, естественно, должно быть равно числу к элементов конструкции. Мат-

рица [Ь] размерностью 4к х т определяет связь между смещениями [д] и дисторсией [Р]. Она со-

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

Полярное разложение дисторсии. Тензор дисторсии представляет скалярное произведение симметричного тензора V и ортогонального R:

F = VR (4)

или наоборот, ортогонального тензора R и симметричного U = RT V R. Тензор R - это тензор жесткого поворота тела на угол p против часовой стрелки в плоскости {e1, e2}:

R = (єіЄі+Є2Є2) cosp + (e2el-ele2) sinp (5)

а V - левый тензор растяжения, который определяется тремя координатами V1, V2 и V3:

V = Vieiei + У2ЄіЄ2 + V2e2ei + V3e2e2.

Он равен сумме единичного тензора I и симметричного тензора деформации є в элементе, предварительно повернутом как жесткое целое на угол p. Тензор U является правым тензором растяжения, он равен сумме единичного тензора и тензора деформации в элементе до его жесткого поворота. В данной работе используется левый тензор растяжения V, с помощью которого проще определить деформации в элементе после его жесткого поворота:

є = V - I = FRT - I. (б)

Такой тензор понадобится позже для расчета напряжений и усилий, приложенных к узлам конструкции.

Выражение (4) представляет двухвалентное тензорное уравнение. Для плоской задачи при известных четырех координатах тензора дисторсии оно содержит четыре неизвестных параметра: p - угол поворота как жесткого целого, однозначно определяющий ортогональный тензор R, согласно выражению (5), а также три координаты симметричного тензора растяжения. В работе [1О] показано, что тангенс, косинус, а также синус угла p вычисляются следующим образом:

tgp-(F3 -F2)/(F1 + F4), cosp-(1+ tg2p)_1/2, sinp-cosp tgp . (7)

Три координаты тензора деформации конечного элемента определяются с использованием выражений (б) и (7):

' Fi'

" Є1" "1" cosp - sinp О О "

Є2 -[Ri] F2 F3 - О , [Ri]- sinp/2 cosp/2 cosp/2 - sinp/2

[Єз _ 1 О О sinp cosp

[F4 _

Матрица [^1] связывает координаты тензора дисторсии элемента с координатами тензора деформации.

Для конструкции тензор деформации записывается, как и тензор дисторсии, в виде столбца, состоящего из k одинаковых блоков:

Ы=[К £2 £з ]1 [^1 £2 £з ]2 .. 1-£1 £2 £з]к Г ,

где е1 и ез - продольные деформации, а 82 - сдвиговая. Выражение, определяющее деформации конструкции, имеет следующий вид:

[8] = Щд)] [Р] - [/2], (8)

где матрица [К(д)]> размерностью зкх4к и состоящая из блоков [^], определяет связь между дис-торсией и деформациями конструкции, а столбец [/2] содержит координаты единичного тензора и имеет размерность зк:

[/2 ]=[[1 0 1] [1 0 1] ... [1 0 1]]Т.

Физическое соотношение модели. В качестве физического соотношения используется закон Гука: о = С-8, где четырехвалентный тензор С является тензором констант упругости. Для конструкции это выражение приобретает матричный вид

[о] = [С] [8], (9)

где столбец [о] размерностью зк содержит компоненты тензоров напряжений оши (о1 и оз -нормальные напряжения, о2 - касательное) элементов конструкции:

[°] = [[°1 0-2 0 ]1 [о 0-2 0 ]2 ... [01 О2 0з ]к ],

а симметричная матрица [С] размерностью ЭкхЭк содержит константы упругости материала. Например, для конечного элемента конструкции из изотропного материала при плоском напряженном состоянии матрица [С] имеет вид:

~ 1

0

[С ]-

E

1 -Az

А

О

1 -А О

где Е - модуль упругости, а # - коэффициент Пуассона.

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

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

рАВ = аАВС • пАВ , рВС = аАВС • пВС, рСА = аАВС • пСА. (10)

Система этих распределенных сил, действующих на элемент, эквивалентна системе сосредоточенных сил 0А, 0В и ЦС, приложенных к узлам корсета (рис. 2, б). Узловые силы 0 определяются с учетом равновесия элемента под действием распределенных сил (АВ рАВ + ВС рВС + СА рСА = 0) и граничных условий (10):

QA --—аABC . nBC, QB --—аABC . nCA,

2 2 '

а) конечный элемент без корсета

б) конечный элемент с корсетом

Рис. 2. Силы, действующие на конечный элемент

Рис. 3. Схема приведения сил к узлу конструкции

Один узел может принадлежать сразу нескольким элементам конструкции, например, узел D, изображенный на рис. 3. Следовательно, сила, приложенная к этому узлу, представляет полусумму сил, действующих на сторонах, не содержащих этого узла. Например, поскольку узел D принадлежит трем конечным элементам DAB, DCA и DBC (рис. 3), сила QD, действующая на него, рассчитывается следующим образом:

Q D = -- (AB * p AB + CA * p CA + BC * p BC) =

= - - (AB * ct DAB ■ n AB + CA * aDCA ■ nCA + BC * aDBC ■ n BC).

2

Для произвольной плоской конструкции это выражение принимает матричный вид:

[Q] = [H(q)] [a], (11)

где столбец [Q] содержит координаты сил, приложенных к узлам конструкции: [Q] = [Q1 Q2 ... Qm]T (по аналогии со смещениями [q]), а матрица [H(q)] связывает координаты тензоров напряжений элементов конструкции с координатами сил, приложенных к узлам конст-

рукции в ее деформированном состоянии, причем эта связь линейна. Матрица [Н(?)] зависит от смещений, так как равновесие элемента рассматривается в деформированном состоянии.

Расчет конструкции при кинематическом нагружении. Результирующее выражение для расчета конструкции следует из решения системы уравнений (3), (8), (9) и (11) относительно нагрузок О]:

О] = №)] [?] + [К?)], (12)

где

№)] = №)] [С] [Вд] [Ь], [%)] = [Н(?)] [С] ( [Вд] [/1] - [/2] ) .

При расчете приложенных к конструкции сил по заданным смещениям история деформирования конструкции не имеет значения.

Расчет конструкции при силовом нагружении. Расчет конструкции при силовом нагружении качественно отличается от расчета при кинематическом нагружении. Во-первых, найти смещения из выражения (12) в явном виде невозможно, так как от смещений зависят матрицы [К(?)] и [£(?)]:

[?] = [К(?)]-1 ( [О] - [Ь(?)] ). (13)

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

Алгоритм расчета

Шаг № 1:

1) задаем силы [О1];

2) принимаем начальное приближение смещений [?нач1] = [0];

3) с помощью уравнения (13) уточняем смещения [?1] на первом шаге.

Шаг № 2:

1) задаем силы О2] = [О1] + [ёО], где [ёО] - приращение сил;

2) задаем начальное приближение смещений [?нач2] = [?1];

3) решаем уравнение (13), уточняя смещения [?2] на втором шаге.

Расчет последующих шагов выполняется аналогично.

Примеры расчетов

Кинематическое нагружение балки при чистом изгибе. В работе выполнены расчеты изгиба балки (рис. 4) при кинематическом нагружении. Этот и все дальнейшие расчеты с применением предложенной расчетной модели выполнены в пакете МЛТЬЛБ. В расчетах использовали следующие исходные данные: материал балки является изотропным и линейно-упругим (модуль упругости Е = 2105 МПа, коэффициент Пуассона # = 0). Длина балки I = 15 м, высота Н = 8 м, толщина 1 м, взаимный угол поворота торцов а = п/2, напряженное состояние является плоским. На рис. 4 приведена схема балки в начальном (с разбивкой на конечные элементы) и в конечном состояниях, число конечных элементов - 240. Координаты узлов сетки вычисляли, исходя из следующих соображений:

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

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

в) длина нейтрального слоя в деформированном состоянии равна начальной длине балки I;

г) поперечные сечения остаются плоскими (гипотеза Бернулли).

Для примера координаты узла А0 в деформированном состоянии балки (точка А на рис. 4) равны:

х = (р - У0) эт^, у = р - (р - Ус) есв^, где угол +, равный хо /р, представляет угол поворота поперечного сечения, содержащего точку А.

В расчете получено распределение сил, приложенных к узлам балки, вполне соответствующее чистому изгибу. На рис. 5 показаны силы, действующие в узлах сетки на правом торце балки, прерывистая линия соответствует известному из курса сопротивления материалов распределению нормальных напряжений по высоте поперечного сечения. На левом торце распределение узловых сил является аналогичным, а в остальных узлах величины внешних сил оказались примерно на два порядка меньшими, чем максимальные, что говорит об относительной корректности решения. Система сил, изображенных на рис. 5, эквивалентна моменту пары сил М, представляющему сумму произведений этих сил на соответствующие расстояния до точки С. В нашем случае моменты на торцах балки составили М = 9,39* 1011 Нм.

В работе выполнено сравнение полученного решения с точным решением, известным из курса сопротивления материалов М = ЕН3а/(12/).

Это выражение соответствует решению задачи о чистом изгибе балки при больших перемещениях в случае равенства нулю коэффициента Пуассона. Вычисленные моменты были сопоставлены с величинами моментов, полученными из расчетов с применением пакета прикладных программ А^УБ. В пакете А^УБ были выполнены два расчета: в первом использованы треугольные конечные элементы (разбивка на элементы совпадает с показанной на рис. 4), во втором - квадратные той же высоты. В расчетах применяли элементы типа р1апе182. Результаты сравнения моментов, вычисленных различными методами, и соответствующих погрешностей , по сравнению с точным решением приведены в таблице. Расчеты показали хорошую адекватность предложенной нами расчетной модели: при выбранном размере элементов погрешность по сравнению с точным решением составляет 5,13 %, что в 4-6 раз меньше, чем в расчетах с применением пакета А^УБ. При уменьшении угла а погрешность предложенного решения остается примерно на том же уровне (около 5 %), а с ростом числа элементов величина , падает. Например, при уменьшении размера элемента вдвое погрешность в расчете с помощью предложенной модели составила 1,2 %, а в расчете с применением пакета А^УБ величина , не изменилась (при увеличении размера элемента вдвое погрешность расчета с применением пакета АК БУБ возросла до 23,6 %).

Таблица

Точное решение Полученное решение Решение с применением пакета ANSYS

Треугольные элементы Прямоугольные элементы

Значение момента, M10-11 Нм 8,94 9,З9 11,4 10,8

Погрешность по сравнению с точным решением, , % - 5,1З 28,0 20,6

Силовое нагружение балки при чистом изгибе. С использованием предложенной модели в работе выполнен расчет конструкции при силовом нагружении. Объект расчета - консольная балка размерами lxh, исходные данные - те же, что и в предыдущем расчете за исключением высоты балки h, которая была уменьшена вдвое. Для решения системы нелинейных уравнений (16) использовали trust region dogleg method, встроенный в функцию fsolve пакета MATLAB. На рис. 6

Рис. б. Узловые силы на правом торце балки

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

V = 4О/3/(ЕН3). Расчет показал заметную нелинейность зависимости прогиба от нагрузки. Например, в случае, когда величина максимального прогиба V оказывается равной высоте балки Н, разница между найденным прогибом и прогибом, вычисленным по линейной теории, составляет около 30 %.

Неоднозначность решения. Геометрическая нелинейность задачи в некоторых случаях приводит к неоднозначности решения: одним и тем же нагрузкам О], приложенным к конструкции, в зависимости от истории изменения нагрузок могут отвечать различные перемещения [?]. Например, для элемента АВС (Е = 2105 МПа, # = 0, толщина 1 м), изображенного на рис. 7, а и б (начальное положение элемента показано сплошной линией), одним и тем же нагрузкам ОАтах = 1,78-1011 Н и Овтах = 2,42-1011 Н отвечают два его различных положения, показанные пунктирными линиями. В первое положение (?А = 4 м, ?в = 0) элемент попадает в случае, когда история изменения нагрузок ОА и ОВ описывается функциями, графики которых изображены на рис. 7, в. Если история нагружения соответствует графикам, изображенным на рис. 7, г, тогда элемент из начального положения попадает в другое положение (?А* = 1,37 м, ?в* = 1,99 м). В первом расчете (рис. 7, в) пропорционально изменяли смещения, а во втором

Рис. 6. Силовое нагружение консольной балки

(рис. 7, г) - силы.

а) элемент в начальном положении и в положении № 1

б) элемент в начальном положении и в положении № 2

в) история изменения нагрузок № 1

г) история изменения нагрузок № 2

Рис. 7. Два варианта истории изменения нагрузок элемента АВС

Заключение

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

Сопоставление с точным решением задачи о чистом изгибе балки при больших перемещениях показывает, что при относительно небольшом количестве конечных элементов ошибка расчета по сравнению с точным решением не превышает 5 %, причем при уменьшении размера элементов ошибка уменьшается. Для сравнения: ошибка расчета в пакете ANSYS при том же размере элементов составляет 28 % для треугольных конечных элементов и около 20 % для квадратных.

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

Литература

1. Трусделл, К. Первоначальный курс рациональной механики сплошной среды / К. Трусделл. - М: Мир, 1975. - 592 c.

2. Chadwick, P. Continuum mechanics: concise theory and problems / P. Chadwick. - 2 ed. -Dover publications, 1998. - 193 р.

3. Belytschko, T. Nonlinear finite elements for continua and structures / T. Belytschko, W.K. Lin, B. Moran. New York: John Wiley and sons, 2000. - 660 p.

4. Мейз, Дж. Теория и задачи механики сплошных сред / Дж. Мейз. - М.: Мир, 1974. - 318 с.

5. An objective time-integration procedure for isotropic rate-independent rate-dependent elastic-

plastic constitutive equations / G.G. Weber, A.M. Lush, T.A. Zavaliangos, L. Anand // International

journal of plasticity. - 1990. - V. б. - P. 701-744.

6. ANSYS - а general purpose finite element program. Rev. 5.0. - Houston (PA): Swanson analysis system inc., 1996. - 510 р.

7. Садаков О.С. онечные деформации в механике деформируемого твердого тела // Вестник ЮУрГУ. Серия «Математика, физика, химия». - 2005. - Вып. б. - 86(46). - С. 114-121.

8. Douglas, R.G. On majorization, factorization, and range inclusion of operators on Hilbert space / R.G. Douglas // Proc. Amer. math. Soc. - 1966. - 8 17. - P. 413-415.

9. Sobczyk, G. Hyperbolic number plane / G. Sobczyk // College mathematics journal. - 1995. -8 26. - P. 268-280.

10. Щербакова, А.О. Использование круга Мора для решения задачи полярной декомпозиции

при плоском напряженном состоянии / А.О. Щербакова, О.С. Садаков, С.И. Шульженко // Вестник ЮУрГУ. Серия «Математика. Механика. Физика». - 2010. - Вып. 2. - 8 9(185). - С. 21-26.

Поступила в редакцию 2 мая 2011 г.

THE FINITE ELEMENT METHOD USE FOR LARGE DISPLACEMENTS CALCULATIONS OF PLANE LINEAR-ELASTIC STRUCTURES

A.O. Scherbakova'

A proposed version of the finite element method can find large displacements of a plane linear-elastic structure with given nodal loads, and conversely, nodal loads at the given displacements. A calculation model includes geometrical, static and physical relations (Hooke's law). The calculation showed that several different deformed positions of the structure may correspond to the same loads. To verify the proposed model, the pure bending of a beam was calculated with large displacements taken into account. Analysis of the results showed the adequacy of the proposed model. At that the results differ markedly from the ANSYS solution.

Keywords: geometric nonlinearity, large displacements, finite deformations, pure bending, plane stress-strain state.

References

1. Truesdell C.A. A first course in rational continuum mechanics. Baltimore, The Johns Hopkins University, 1972.

2. Chadwick P. Continuum mechanics: concise theory and problems (2nd ed.). Dover publications, 1998. 193 p.

3. Belytschko T., Lin W.K., Moran B. Nonlinear finite elements for continua and structures. New York, John Wiley and sons, 2000. ISBN 0-471-98773-5. 660 p.

4. Mase G.E. Theory and problems of continuum mechanics: Schaum's Outline Series. New-York: Mcgraw-hill Book Company, 1970. 221 p.

5. Weber G.G., Lush A.M., Zavaliangos T.A., Anand L. An objective time-integration procedure for isotropic rate-independent rate-dependent elastic-plastic constitutive equations. International journal ofplasticity. 1990. Vol. 6. pp. 701-744.

6. ANSYS - a general purpose finite element program. Rev. 5.0. Houston (PA): Swanson analysis system inc., 1996. 510 p.

7. Sadakov O.S. Vestnik YuUrGU, serija «Matematika, fizika, khimija». 2005. Vol. 6, no. 6(46). pp. 114-121. (in Russ.).

8. Douglas R.G. On majorization, factorization, and range inclusion of operators on Hilbert space. Proc. Amer. math. Soc. 1966. no 17. pp. 413-415.

9. Sobczyk G. Hyperbolic number plane. College mathematics journal. 1995. no. 26. pp. 268-280.

10. Scherbakova A.O., Sadakov O.S., Shul'zhenko S.I. Ispol'zovanie kruga Mora dlja reshenija zad-

achi poljarnoj dekompozicii pri ploskom naprjazhennom sostojanii (Application of the Mohr circles for solving the problem of polar decomposition under plane stress). Vestnik YuUrGU. Serija «Matematika. Mehanika. Fizika». 2010. Vol. 2, no. 9(185). pp. 21-26. (in Russ.).

1 Scherbakova Alla Olegovna is Cand. Sc. (Engineering), Applied Mechanics, Dynamics and Strength of Machines Department, South Ural State^niversitye-mail:^llaScherbakova@list.ru^^^^^^^^^^^^^^^^^^^^^^^^^_^^^^^^^^^^^^^^^^^^^^^^^^^_