21. Хохлов А.В. Свойства кривых ползучести при ступенчатом пагружепии линейного определяющего соотношения вязкоупругости // Пробл. прочности и пластичности. 2015. 77, № 4. 344-359.
22. Шестериков С.А., Юмашева М.А. Конкретизация уравнения состояния при ползучести // Изв. АН СССР. Механ. твердого тела. 1984. № 1. 86-91.
Поступила в редакцию 25.05.2016
УДК 539.3
О ТЕНЗОРНЫХ МЕРАХ НАПРЯЖЕНИЙ И ДЕФОРМАЦИЙ, ИСПОЛЬЗУЕМЫХ В ANSYS ДЛЯ РЕШЕНИЯ УПРУГОПЛАСТИЧЕСКИХ ЗАДАЧ ПРИ КОНЕЧНЫХ ДЕФОРМАЦИЯХ
Н. В. Овчинникова1
Установлено, что при решении упругопластических задач в пакете ANSYS в качестве объективной производной для тензора напряжений Коши использована производная Динса (Грина-Мак-Инниса-Нахди), а не производная Яуманна, как заявлено в теоретическом руководстве к ANSYS. Вычисляемый для этих задач тензор деформаций является не тензором логарифмических деформаций Генки, как указано в теоретическом руководстве к ANSYS, а правым неголономным тензором деформаций, таким, что связанный с ним левый неголономный тензор деформаций порожден производной Динса.
Ключевые слова: упругопластическое тело, конечные деформации, меры деформаций, меры напряжений, пакет ANSYS.
It is shown that the Dienes (or the Green-Mclnnis-Naghdi) stress rate is used as an objective stress rate of the Cauchy stress for solving elastic-plastic problems in the ANSYS program instead of the Jaumann stress rate announced in the ANSYS Theory Reference. The strain tensor is not the logarithmic or Hencky strain measure as is announced in the ANSYS Theory Reference, but it is a right nonholonomic strain tensor such that the left tensor is generated by the Dienes derivative.
Key words: elastoplastic body, finite strain, strain measures, stress measures, ANSYS program.
В теоретическом руководстве fl] к конечно-элементному пакету ANSYS указано, что вычисление тензора деформаций при решении задач для упругопластических тел в случае конечных деформаций производится с помощью алгоритма, предложенного Т. Хьюзом в работе [2]. В результате сравнения описания алгоритма вычисления деформаций, приведенного в [1], с представленным в работе [2], а также аналитического и численного решения в ANSYS задачи о простом сдвиге гипо-упругой среды нами установлено несоответствие заявленных в [1] и используемых в расчетах мер напряжений и деформаций.
Пусть движение среды описывается законом х% = ж^, t), где t — время; Xi, Xj (i,j = 1,2,3) — координаты (в декартовой системе координат) точки в актуальной и отсчетной конфигурациях. Пусть F — аффинор деформации: Fij = dxi/dxJ = |detF|. По теореме о полярном разложении F = RU = VR, где R — ортогональный тензор (тензор поворота), U — правый и V — левый тензоры искажений.
Приведем описание алгоритма пошагового вычисления деформаций в соответствии с документацией к ANSYS [1]. Пусть для всех шагов по времени вплоть до п-го шага tn известны координаты точки тела, т.е. известны векторы х°,...,хп, где xk = x(x°,tk). Следовательно, известны векторы и0,... ,ип перемещений точки: ик = хк — х°. Пусть для всех шагов вплоть до (п — 1)-го шага вычислены значения компонент тензора "деформаций" е. Тогда для того чтобы определить значения компонент тензора "деформаций" еп для п-го шага, требуется:
1 Овчинникова Нелли Викторовна — канд. физ.-мат. наук, доцент каф. теории упругости мех.-мат. ф-та МГУ, e-mail: ovch-nQyandex.ru.
1) вычислить "среднее значение перемещений" точки для (п — 1)-го и п-го шагов:
и0'5 = (и"-1 +ип)/2; (1)
2) вычислить радиус-вектор точки, соответствующий этому "среднему перемещению":
х0'5 =х°+ и0'5; (2)
3) вычислить приращение перемещений:
А и11 = (и11 - и11-1)] (3)
4) вычислить тензор Ае" с компонентами
1 ( дА и? дАи
Ае-- =
5) определить аффинор соответствующий "среднему перемещению":
(4)
+ ь-? = 1,2,3; (5)
6) найти полярное разложение аффинора
Г""' /Г'"-(6)
7) вычислить приращение деформаций Аеп:
Аеп = (Е°'5)Т • Аеп • (Е0'5) ; (7)
8) вычислить полное значение деформации еп для п-го шага:
еп = еп~1 + Аеп. (8)
Алгоритм (1)-(8) из теоретического руководства к А^УЭ отражает представленный в работе [2] алгоритм вычисления деформаций. Но в руководстве к А^УЭ утверждается, что по данному алгоритму вычисляются компоненты тензора логарифмических деформаций Генки О = 1п [7, а Т. Хьюз в своей работе рассматривал способ вычисления другой меры деформаций — некоторого тензора, обозначенного в [2, ф-ла (81)] через ей, такого, что
ен = КТУК, (9)
где ¿я — производная по времени тензора ей, V = (с1 + с1т)/2 — тензор скоростей деформации, с1 — тензор скоростей дисторсий с компонентами = дгц/дх^ {иг = дХг/&Ь — компоненты вектора и скорости движения точки). Определение (9) эквивалентно следующему (пусть вд(£о) = £о — начальное значение тензора деформаций):
£д(£)=£о + J КТУКйт. (10)
¿0
Сопоставляя (10) с (1)-(8), видим, что тензор еп (8) действительно представляет собой приближение к тензору £д(£га).
Тензор отвечающий определению (9) или (10), известен в литературе [3-5] как правый (материально ориентированный) неголономный тензор деформации ("правая скоростная мера деформации" [3]). Согласно [3, 5], связанный с ним левый (пространственно ориентированный) тензор ег.
ег = КекКт (11)
является неголономным тензором конечных деформаций, порожденным производной Динса ("левой скоростной мерой деформации" [3]), т.е. в каждый момент времени тензор £1 является решением уравнения
В{ег) = V. (12)
Здесь 0{е{) — производная Динса [6] тензора £1 (которая также известна в литературе как производная Нахди-Грина-Динса [2], Грина-Мак-Инниса-Нахди [7], "нейтральная" производная [3]), определяемая следующим образом:
В(£1) = (¿! - П0е1 + £1П о), (13)
где ¿1 — полная (материальная) производная тензора £1 по времени, Г^о — ИИ •
Следует отметить, что левый неголономный тензор конечных деформаций, связанный с £д (9) соотношением (11), вводится в работе [2, ф-ла (52)] как тензор £1, отвечающий уравнению ¿1 = V, что является ошибкой, поскольку £1 должен удовлетворять уравнению (12). Но этот тензор носит в [2] вспомогательный характер и не участвует в алгоритме вычисления напряжений.
Известно [3], что правая ед и левая £1 неголономные меры деформаций совпадают соответственно с правым С = 1п и и левым II = 1п V тензорами логарифмических деформаций Генки только для тех движений, когда главные оси деформации (тензора II) фиксированы во времени. Для произвольных же случаев движения эти тензоры различны, т.е. £д ф О, £1 ф Н. Примером движения, при котором главные оси тензора деформации в каждый момент времени определяются разными материальными волокнами, является простой сдвиг среды. Тензор £д для этого движения может существенно отличаться от О [8, 9].
Рассмотрим простой сдвиг среды, когда движение задано в декартовой системе координат уравнениями х\ = х®, Х2 = х^ + к-х®, жз = Ж3) где к = ко1 — величина сдвига. Примем величину скорости сдвига ко равной 1, т.е. к =
Значения ненулевых компонент правого тензора деформаций Генки имеют вид
,___Ф-— (14)
Вычисленные аналитически по алгоритму (1)-(8) значения ненулевых компонент тензора деформаций £п при заданном количестве шагов расчета п для достижения заданной величины сдвига I (шаг расчета соответственно равен = 1/п) определяются выражениями
т-1 (15)
1Д 4 — (т — 0,5)2А^2 д
т= 1
Значения деформаций (15) совпадают со значениями компонент тензора деформаций, полученными при решении в А^УЭ задачи о простом сдвиге с использованием элементов Р1.Л.\Т. 12. Р1.Л.\Т.1К2. Р1.Л.\К1 КЗ при заданных ¿и п. Это подтверждает тот факт, что в качестве алгоритма вычисления конечных деформаций в А^УЭ принят алгоритм (1)-(8). Заменяя в (15) конечную сумму интегралом, получим предельные при —> 0 значения компонент е^, равные компонентам тензора £д(£):
I
= ~ек33(1) = -2 У йт = -1п(1 +¿2/4),
(16)
1 [А -г2 = 2 у йт = + 2аг<*е(*/2).
о
Сделав проверку, можно убедиться, что тензор ед с компонентами (16) действительно является решением уравнения ¿д = КТУЯ, т.е. отвечает определению (9).
Сравнение значений компонент тензора деформаций (16) и тензора логарифмических деформаций Генки (14) при 0 ^ I ^ 2тт представлено на рисунке. Видно, что компоненты тензора еRij близки к деформациям Генки О^ только при I < 1 (когда \ещ\ < 0,4). С ростом I отличие компонент тензоров ещ и О^ быстро растет и становится существенным.
Сравнение; компонент Goo и £r22 (о). G03 и ед23 (б) Зная компоненты тензора еr и компоненты тензора поворота II
Rii — 1/Vt2 + 4, R22 — R33 — 2-Rn, R23 — —R32 — tRu, можно вычислить компоненты тензора t/ = RerRt:
t/22 t-i зз
4 - t2 , 4 + i
In-
2i2
4i , 4 + i2 In —■—
4 + i2 " 4
i 4-i2
+
8i
4 + i2 1 4 + i2arCtg2'
4 + i2
4 — i2 i
- -0 + 2 --гт arct.g -.
24 +i2 4 + i2 2
(17)
Сделав проверку, можно убедиться, что тензор е/ с компонентами (17) действительно является решением уравнения D{ei) = V.
Таким образом, на примере задачи о простом сдвиге проверено, что вычисляемый в ANSYS тензор деформаций является правой неголономной мерой деформации £r, а связанный с ним левый тензор t/ порожден производной Динса. Зная это, легко установить, что в качестве объективной производной для тензора напряжений Коши в принятых в ANSYS определяющих соотношениях теории пластического течения при конечных деформациях используется не производная Яуманна (как утверждается в [1|), а производная Динса. Для этого достаточно решить в ANSYS задачу о простом сдвиге для гииоуиругой модели. Действительно, пусть в ANSYS в качестве определяющего соотношения теории пластического течения при конечных деформациях принято соотношение, предложенное в работе [2|:
D{a) = a:{V-Vp), (18)
где D(a) производная Динса (13) для тензора напряжений Коши а; а : V свертка тензора а (4-го ранга) упругих постоянных материала с тензором (V — Vp), где V тензор скоростей деформаций; Vp пластическая часть тензора скоростей деформации. Определяющее соотношение (18), записанное для левых тензоров а и V, эквивалентно [2| соотношению для правых тензоров аr и RTVR:
&r = А : (RTVR - RTVVR) , (19)
где Gr = RTaR "повернутый" [7] тензор напряжений Коши, дr производная по времени тензора ctr-
Пусть пластическая часть тензора скоростей деформации равна нулю, т.е. в (18) и (19) имеем Vp = 0. Рассмотрим случай изотропной среды, когда компоненты тензора а в декартовой системе координат имеют вид а^д./ = Лöijöki + M^fc^j/ + öuöjk), где Л, /л константы Ламе, öjk символ
Кроиекера. Тогда соотношение (18) или (19) становится соотношением гипоупругости вида D{a) = а : V или
àR = a: (RTVR). (20)
С учетом (9) соотношение (20) можно переписать следующим образом:
&r = a: £R. (21)
Проинтегрировав (21) по времени, получим соотношение, связывающее правые неголономные меры напряжений и деформаций: <jr = а : er или ar = А ■ tr{ев)I + 2¡J£r, где I — единичный тензор, tr{er) — след тензора er.
Для задачи о простом сдвиге ít{er) = 0, следовательно, a r = 2 ¡jler и компоненты тензора a r имеют вид
(TR22 (í) = -<TR33 (í) = -2/x ln(l + t2/4), CFR23{t) = 2/x(—0,5í + 2arctg(í/2)).
Поскольку a = RorR , то a = 2 ¡iRerRt = 2 ¡l£h следовательно, компоненты тензора а имеют вид
/ 4 — i2, 4 + í2 2 t2 8t t'
= -assit) = 2/x In ■— - ^ + ^ arctg
( 4í , 4 + í2 t 4 — t2 4 — t2 t4 (23)
"23(f) = 2/x ^^ ln —4 2 4 + í2 + 4 + i2 arCtg 2,
Можно показать, что тензор напряжений а с компонентами (23) совпадает с решением задачи о простом сдвиге изотропной гипоупругой среды, т.е. с решением уравнения D{a) = 2fxV, полученным Дж. Динсом в работе [6].
Определяющее соотношение для гипоупругой среды с использованием производной Яуманна записывается следующим образом:
DJ{a) = а : V, (24)
где DJ{a) = à — Wa + aW — производная Яуманна (производная Зарембы-Яуманна-Нолла) для тензора a, W = {d — dT)/2 — тензор скоростей вращений. Решение задачи о простом сдвиге гипоупругой среды для определяющего соотношения (24) имеет вид [6, 10]
022 (t) = -0-33 (t) = ju(l — cosí), ^
0"23 {t) = Ai sin í.
Напряжения (23), полученные для определяющего соотношения с использованием производной Динса, и напряжения (25), полученные для определяющего соотношения с использованием производной Яуманна, близки только для величины сдвига t < 1, а далее они существенно различаются. В частности, напряжение 023{t), определенное формулой (23), является монотонно возрастающей функцией, а напряжение 023{t), определенное формулой (25), — гармонической функцией.
Решение в ANSYS задачи о простом сдвиге для упругого материала в случае величины сдвига 0 < t ^ 2тт показывает, что вычисленные в ANSYS значения компонент тензора напряжений Коши совпадают с <Tíj{t), выраженными формулами (23). Графики зависимости компонент "повернутого" тензора напряжений от величины сдвига t совпадают с графиками компонент тензора o\r(í), построенными по формулам (22). (Эти графики с точностью до коэффициента 2¡jl совпадают с графиками компонент er, представленными на рисунке.) Следовательно, при решении упругопластических задач в ANSYS используются определяющие соотношения именно вида (18), т.е. с применением производной Динса, а не производной Яуманна.
Замечание. Для большинства элементов в ANSYS доступен вывод значений компонент как левых (£i,cr), так и правых {er,ctr) тензоров напряжений и деформаций. Компоненты тензора напряжений Коши а относительно базиса {e¿, г = 1,2,3} декартовой системы координат совпадают с компонентами тензора a r относительно "повернутого" базиса {üe¿}, т.е. выполнено а = a^eiej, cr = (тíjReíRej. Для представления в постпроцессоре (General Postproc) значений компонент относительно базиса e¿ правых тензоров деформаций er и "повернутого" тензора напряжений or требуется предварительно выполнить команду
Main menu: General Postproc —>• Options for Output —>• RSYS = As calculated.
Для вывода значений компонент относительно базиса {е^} левых тензоров деформаций ei и напряжения Коши а требуется предварительно выполнить команду
Main menu: General Postproc —> Options for Output —> RSYS = Global Cartesian.
В TimeHistory Postproc (постпроцессоре истории нагружения) строятся графики зависимости от "времени" только для правых тензоров er и <jr.
Работа выполнена при поддержке грантов РФФИ № 15-08-04281, 16-01-00669.
СПИСОК ЛИТЕРАТУРЫ
1. ANSYS Academic Teaching Introductory, Release 16.2, Help System, Mechnical APDL, Theory Reference. Ch. 3. Structures with Geometric Nonlinearities, ANSYS, Inc., 2015.
2. Hughes T.J.R. Numerical implementation of constitutive models: rate-independent deviatoric plasticity // Theoretical foundation for large-scale computations for nonlinear material behavior / Ed. by S. Nemat-Nasser, R.J. Asaro, G.A. Hegemier. Dordrecht, Netherlands: Martinus Nijhoff Publishers, 1984.
3. Бровко Г. JI. Некоторые подходы к построению определяющих соотношений пластичности при больших деформациях // Упругость и неупругость. М.: Изд-во МГУ, 1987. 68-81.
4. Толоконников Л. А., Маркин А.А. Определяющие соотношения при конечных деформациях // Проблемы механики деформируемого твердого тела: Межвуз. сб. тр. Калинин, политех, ин-т. Калинин: Изд-во КГУ, 1986. 49-57
5. Маркин А.А., Соколова М.Ю. Термомеханика упругопластического деформирования. М.: Физматлит, 2013.
6. Dienes J.К. On the analysis of rotation and stress rate in deforming bodies // Acta Mech. 1979. 32. 217-232.
7. Simo J.C., Hughes T.J.R. Computational Inelasticity. NY: Springer-Verlag, 1998.
8. Овчинникова H.B., Муравлев А.В. Оценка точности вычисления инкрементальными методами конечных упругопластических деформаций в пакетах ANSYS и QFORM7 //XI Всерос. съезд по фундаментальным проблемам теоретической и прикладной механики: Сб. тр. (Казань, 20-24 августа 2015 г.). Казань: Изд-во Казан. (Приволжского) федер. ун-та, 2015. 2825-2827.
9. Овчинникова Н.В., Муравлев А.В. О применении мер конечных деформаций в пакетах ANSYS и QFORM7 // Упругость и неупругость: Мат-лы Междунар. науч. симп. по проблемам механики деформируемых тел, посвященного 105-летию со дня рождения А.А. Ильюшина (г. Москва, 20-21 января 2016). М.: Изд-во МГУ, 2016. 233-236.
10. Johnson G.C., Bammann D.J. A discussion of stress rates in finite deformation problems // Int. J. Solids and Struct. 1984. 20, N 8. 725-737.
Поступила в редакцию 12.09.2016
УДК 539.3
СИНТЕЗ НЕОДНОРОДНОЙ УПРУГОЙ СИСТЕМЫ С ГРАНИЧНОЙ НАГРУЗКОЙ
Л.Д. Акуленко1, A.A. Гавриков2, C.B. Нестеров3
Для адекватного описания свободных и управляемых движений одномерных упругих систем с распределенными параметрами рассматривается соответствующая модель, выражаемая линейной краевой задачей с граничными условиями типа условий третьего рода. Считается, что управляющее воздействие входит аддитивно в уравнение движения
1 Акуленко Леонид Денисович — доктор физ.-мат. наук, гл. науч. сотр. ИПМех им. А.Ю. Ишлинского РАН, e-mail: l.akulenko®bk.ru.
2 Гавриков Александр Александрович — канд. физ.-мат. наук, ст. науч. сотр. ИПМех им. А.Ю. Ишлинского РАН, e-mail: gavrikovQipmnet.ru.
3 Нестеров Сергей Владимирович — доктор физ.-мат. наук, проф., гл. науч. сотр. ИПМех им. А.Ю. Ишлинского РАН, e-mail: kumakQipmnet.ru.