CC BY
40
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
Том 14 7, ки. 3
Физико-математические пауки
2005
УДК 539.3
большие вязкоупругопластические деформации трехмерных тел
А.И. Голованов, Л. У. Султанов
Аннотация
Предложен метод исследования напряженно-деформированного состояния (НДС) вяз-коупругопластических тел с учетом больших перемещений, поворотов и конечных деформаций. Используется процедура пошагового пагружепия в рамках комбинированного лаграпжево-эйлерового описания деформирования среды. Решение находится с использованием метода конечных элементов (МКЭ).
1. Разрешающее уравнение. Определяющие соотношения
Традиционно в механике деформируемого твердого тела (МДТТ) для решения нелинейных задач получило распространение лагранжево описание среды, согласно которому состояние элементарного объема описывается в компонентах вектора перемещений из недеформированного состояния в деформированное и второго тензора напряжений Пиолы Кирхгофа, также отнесенного к недеформированно-му объему. В этом случае хорошо формулируется краевая задача в дифференциальной или вариационной форме, для решения которой возможно использование различных численных методов. Однако подобный подход имеет существенный недостаток в задачах с большими деформациями. Он связан со сложностью построения определяющих соотношений между используемыми тензорами напряжений и деформаций, которая особенно сильно проявляется при постулировании определяющих соотношений в дифференциальной (скоростной) форме. При этом, если течение среды описывать в эйлеровой постановке, то можно эти трудности обойти. Для описания процесса деформирования используется комбинированная лагранжево-эйлерова постановка. Поведение материальной точки (элементарного объема) отслеживается в соответствии с лагранжевым методом описания среды, а в текущий момент времени процесс деформирования представляет собой течение среды с некоторыми физико-механическими свойствами. Это соответствует подходу Эйлера, широко применяемому в механике сплошных сред.
Для описания процесса деформирования среды вводятся радиус-векторы материальной частицы И. = угег в деформированном состоянии иг = хг ег в неде-формированном состоянии, вектор и = И — г перемещения точки тела из недеформированного в деформированное состояние, вектор V = И = угег скорости
где ег - орты глобальной декартовой системы координат. Напряженное состояние описывается тензором истинных напряжений Коши-Эйлера (Я) = а3ег < е3- . Исходным является уравнение виртуальных мощностей в актуальной конфигурации
перемещения точки, тензор деформаций скорости
JJJ (£) • • (Sd) dQ
III
f • Sv dQ +
II
p • Sv dS,
(1)
где П - текущий объем, - часть его поверхности, на которой заданы усилия, f, р - плотности векторов массовых и поверхностных сил. Кинематические граничные условия ¿V = 0 на поверхности, не принадлежащей П, выполняются априори. Путем линеаризации уравнения мощностей (1) получено уравнение в скоростях напряжений
С •• (6() + (Е)
61
+ >
■(61)
(П
• ¿V (1П
.1 Р1
• ¿V (Б, (2)
.
Помимо скорости изменения напряжений (1С) пепользуется индифферентная производная тензора напряжений по Яуманну Коши- Эйлера (С7), которая имеет
(е7) = (ё) - и • (£) + (£) • и, (3)
1 / д? V \
где = - I —- — \ е?. ej тензор скоростей поворота. Рассматривается
класс нелинейно упругих материалов, для которых предполагается существование потенциальной энергии упругой деформации Ш как функции инвариантов тех или иных мер деформаций, производная от которой по соответствующей мере деформации определяет тензор напряжений, в частности,
/сЖЛ
' \дв)
(4)
где (В) - мера деформации Фингера (левый тензор Коши-Грина). В качестве примера приведен стандартный материал второго порядка [4], для которого определяющее соотношение в скоростях напряжений записывается в следующем виде:
£ =
i *
] I 2
(I) •• В • (В) +
- 3) - м
В +2Д(В) • (В
-^(Е)- (5)
Здесь А, А - коэффициенты Ляме обычного закона Гука, 1\в 12 в - первый и
(В)
При решении задач с учетом пластического деформирования используется аддитивное представление для деформации скорости, считается, что относительное изменение объема является упругой деформацией скорости, предполагается справедливость ассоциированного закона течения. Рассмотрен идеально пластический изотропный материал. Для удобства тензоры деформации скорости и напряжений
разложены на шаровые части ¿о = а о = -а" и девиаторы = № — 6.^(1, о,
а'гз = агз — 6^ а0, где 6ц - символ Кронекера.
Физические соотношения упругого деформирования записаны в виде линейной зависимости между производной Яуманна тензора напряжения (3) и тензора деформации скорости
(Е£) =3К ((о) , (С7) =2С (('),
(6)
где О - модуль сдвига, К - модуль объемного сжатия. В этом случае соотношения упругости будут полностью удовлетворять принципу индифферентности. В качестве условия пластичности в работе используется критерий Губера Мизеса, для которого функция текучести Ф имеет вид
Ф (а3) = аг — ат = 0, (7)
где аг - интенсивность напряжений, ат — передел текучести. Для материалов с внутренним трением применяется условие пластичности Мизеса Боткина
Ф (а«) = - С* + о"о ^ V* =0. (8)
Здесь С % ф* - сцепление и угол внутреннего трения на октаэдрических плогцад-
В рамках классической теории течения пластическое деформирование моделируется на основе метода проецирования напряжений на поверхность текучести [3], т. е., если имеется два бесконечно близких состояния /и (/ +1), то по известным параметрам 1-го состояния определяются напряжения (/ + 1)-го состояния. Сначала
вычисляется тензор пробных напряжений ^Ё^ = ('Я) + ('Ё^ где Ё^ определяется с помощью (3) и (6). Тогда шаровая часть тензора истинных напряжений будет иметь вид
('+%) = (Ёо) , О)
а девиатор истинных напряжений определяется как проекция девиатора «пробных» на поверхность текучести
С+1Е') = £1 (£>\ . (Ю)
аг ^ ^
При моделировании процессов деформирования с учетом деформаций ползучести используется аддитивное представление деформации скорости. Рассмотрена обобщенная модель Максвелла, при учете пластических деформаций используются теория течения и метод проецирования напряжений на поверхность текучести. Определяющие соотношения в этом случае записаны как
(Я') О (и)
Ёо) =3К (¿о),
(Ё') =2О (¿') + М • (Я) — (Я) • И —
где п _ коэффициент вязкости.
Для каждой модели поведения среды получены разрешающие уравнения путем подстановки соответствующих физических соотношений в уравнение (2).
2. Алгоритм решения
Для решения нелинейных задач используется метод последовательных нагру-женпй, который может быть естественно реализован в рамках МКЭ. Процесс деформирования представляется в виде последовательности равновесных состояний. Переход от предыдущего состояния к последующему происходит путем приращения нагрузки. Методика расчета состоит в разработке алгоритма вычисления (/ + 1) /
дом шаге нагруженпя строится разрешающее уравнение, в которое входит особое слагаемое невязка. По физическому смыслу это есть уравнение виртуальных
мощностей (1) в /-м состоянии, которое должно удовлетворяться для точного решения. Так как это уравнение на шаге нагружения в настоящей методике точно не выполняется, то его вводим в линеаризованное уравнение в виде невязки. Опыт решения нелинейных задач в шаговой постановке (в том число и представленной в настоящей работе) свидетельствует, что наличие такого рода слагаемых в правой части соответствующего линеаризованного уравнения препятствует накоплению ошибок н не позволяет решению удаляться от действительной кривой деформирования. Полученное уравнение является линейным относительно скорости 1\. Так как исследуемые процессы не имеют явного динамического характера (ускорения не учитываются), то под временем можно понимать любой монотонно возрастающий параметр, определяющий изменение нагрузки. В таком аспекте вполне уместно принять производную по времени как отношение приращения соответствующих величин, получаемых при переходе с 1-го состояния в (1 + 1)-е, к приращению по времени. В результате конечноэлементной дискретизации получено уравнение для приращений перемещений в узловых точках {Д'м}
[1К] {Д1и[ = {ДгРи - {1Ни , (12)
да [гК] - матрица левых частей, {ДгР} вектор приращения узловых сил, {гЯ} - вектор невязки.
Решая систему алгебраических уравнений (12), получим приращения перемещений, с помощью которых определим конфигурацию (/ + 1)-го шага г+1Н = гН+Дги и напряженное состояние (г+1 Е) = (г Е) + Е^ Дг.
При появлении пластических деформаций и вследствие использования метода проецирования напряжений на поверхность текучести (9), (10) полученное напряженное состояние не удовлетворяет разрешающей системе уравнений (12). Поэтому применяется итерационное уточнение текущего НДС. Эта итерационная процедура основана на введении в разрешающее уравнение вариации мощности; «дополнительных напряжений» (Ед) па возможных деформациях скорости И! С Е™) • • (бг() (П, где дополнительные напряжения определяются как разность довиаторов истинных и «пробных» напряжений. Итоговое уравнение для ш-й итерации па /-м шаге нагружения имеет вид
[гК] {д1пт[ = {ДгРи - с1ни - СгБтИ, (13)
где {гБт} - вектор дополнительных напряжений.
При исследовании вязкоупругопластичоких деформаций алгоритм итерационного уточнения выглядит аналогично, за исключением разрешающего
««« с )
/// % С^') ' ' характеризующее вязкие деформации.
П
При исследовании закритичоского поведения при потере устойчивости нагру-жонио не носит монотонно возрастающего характера, вследствие чего применение классического шагового метода представляется невозможным, тогда используется алгоритм, основанный на методе продолжения по параметру [5]. Процесс деформирования представляется в виде последовательности равновесных состояний, тогда на очередном шаге имеем систему линейных алгебраических уравнений (13). Счи-
/
прикладываомой нагрузки является неизвестной и представляется в следующем виде:
{ДгР} = Д1 г {Р о} ,
где {Р0} - вектор, задающий направление нагружения, а Дгt - искомый параметр, определяющий приращение нагрузки. Тогда решение ищется в виде проекции касательной к кривой деформирования
{Д1п) = 1А {Дг-1и} + {ДУ} , Д^ = 1\Д1-Ч — Дг,
где 1А - параметр, определяющий величину шага, {Дг-1п, Дг-1£} - вектор касательной, {ДУ, Дг} - вектор проекции касательной к кривой деформирования, который находится из условия ортогональности вектора касательной к вектору проекции и уравнения равновесия.
Параметр 'А, определяющий длину касательной на каждом шаге, в приведенных ниже численных расчетах принят как отношение длины дуги кривой деформирования на предыдущем шаге приращения к длине дуги на первом шаге.
Пространственная дискретизация основана на методе конечных элементов в рамках полилинейной трехмерной пзопараметрпческой аппроксимации на базе восьмпузлового элемента.
3. Численные примеры
1. Задача об изгибе в кольцо упругой балки прямоугольного поперечного сечения, жестко защемленной, с одной стороны, и нагруженной изгибающим моментом, с другой. Определяющие соотношения записаны в виде (6). Значение момента, при котором балка изгибается в кольцо, было найдено аналитическим путем. Параметрическое исследование по изменению сотки конечных элементов показало, что увеличение числа элементов по высоте мало сказывается на точности решения, наибольшее влияние оказывает число элементов по длине. Величина шага нагружения существенно влияет на точность, что вполне объяснимо. На рис. 1 изображено деформированное состояние балки и несколько промежуточных этапов нагружения при сетке конечных элементов 200 х 1 х 1 при разбиении нагрузки на 1000 шагов.
Рис. 1. Деформированные состояния балки
2. Распределение напряжений в толстостенной длинной трубе, находящейся под действием осесимметричпого внутреннего давления р при упругопластичпом деформировании в геометрически линейной постановке (плоская задача). Внутренний и внешний радиусы трубы равны а = 1 см и Ь = 2 см соответственно, модуль упругости Е = 2.0-106 кг/см2 , коэффициент Пуассона / = 0.3. Материал полагаем идеально пластическим, критерием пластичности служит условие Губера Мизоса
а/аг
Гис. 2. Гаспредслепне радиальных и окружных напряжений в трубе
(7). физические соотношения упругого деформирования записываются в виде (6). Из аналитического решения было найдено отношение внутреннего давления к пределу текучести р/ат = 0.7208, при котором радиус пластической зоны с =1.5 см. На рис. 2 показано распределение радиальных и окружных напряжений в трубе по отношению к пределу текучести при сетке конечных элементов 80 х 20, а также аналитическое решение (штриховая линия).
Гис. 3. Деформированное состояние плиты
3. Задача об упругом деформировании плиты из стандартного материала второго порядка (5) под действием равномерного давления. Нижнее ребро плиты не имеет вертикального смещения. Плита — квадратная со стороной а = 20 см и толщиной к = 0.5 см, Е = 1000 кг/см2, / = 0.3. На рис. 3 изображено деформированное состояние плиты.
4. Упругопластичоскоо деформирование жестко защемленной с обоих концов балки прямоугольного поперечного сечения под действием распределенной на-
Рис. 4. Распределение интенсивности пластических деформаций в нагруженном и разгруженном состояниях
д, кг/см2
Рис. 5. Зависимость интенсивности напряжений от нагрузки
грузки ц = 25 кг/с м2. Длина бал ки 1 = 25 см, высот а -К =1 см, ширина -Ь = 0.125 »1, Е = 2.0 • 104 кг/см2, / = 0, предел текучести <гт = 750 кг/см2. Материал идеально пластический, подчиняющийся критерию пластичности Губера Мизоса (7), определяющие соотношения заданы в виде (6). Так как задача является симметричной, то достаточно рассмотреть половину балки, введя дополнительные условия симметрии. При решении используем сетку конечных элементов размером 100 х 10 х 1. Нагрузка разбита па 100 шагов.
На рис. 4 показано распределение интенсивности пластических деформаций в нагруженном и разгруженном состояниях для половины балки. Интересно отметить в этой задаче наличие впадины в окрестности защемления, а именно в этой зоне возникают конечные пластические деформации. Наличие этой впадины объясняется учетом физической нелинейности. На рис. 5 представлены зависимости
Н', см
Гис. 6. Зависимости нагрузки от прогиба средней точки арки
интенсивности напряжений от нагрузки в точках, указанных на рис. 4, где сплошной линией показан этап нагружония. а штриховой этап разгрузки. Так. в точке 1 можно наблюдать этап упругого деформирования с возрастанием интенсивности напряжений, затем наступает этап пластического с возрастанием интенсивности пластических деформаций: при разгрузке сначала происходит убывание интенсивности напряжений с неизменными пластическими деформациями, затем увеличение до предела текучести и наступление этапа догрузки с возрастанием пластических деформаций (здесь имеет место влияния впадины). В точках 2 и 4. как видно из графиков, после достижения предела текучести, возникают пластические деформации, которые сохраняются и после снятия нагрузки. Однако в точке 2 можно наблюдать эффект падения интенсивности напряжений при возрастании нагрузки. В точке 3 возникают лишь упругие деформации.
5. Задача устойчивости шарнирно опертой круговой нерастяжимой арки под равномерным давлением. Радиус арки 10 см, центральный угол равен 90° , толщина и ширина арки - 1 см, Е = 2000 кг/см2 , ^ = 0. Физические соотношения выбраны в виде (6).
При решении были получены как симметричные, так и несимметричные относительно середины арки формы деформирования. При симметричной форме рассматривалась половина арки, разбитая на сетку конечных элементов размером 120 х 2 х 1, с постановкой граничных условий шарнирного опирания для конца арки (узлы, расположенные на средней линии основания, полностью защемлены) н условий симметрии для центральной точки (узлы не имеют горизонтального смещения). При несимметричной форме деформирования была рассмотрена целая арка с сеткой конечных элементов 240 х 2 х 1с постановкой граничных условий шарнирного опирания для концов арки. Для того чтобы перейти к несимметричной форме закритичоского деформирования, вводились возмущения на нагрузку. В результате были получены решения для соответствующих форм изгибов, которые дают хорошее совпадение с [5]. Значения критических нагрузок для симметричной д* = 1.343 кг/см2 и несимметричной д* = 0.5 кг/см2 форм практически совпадают с аналитическими решениями дт = 1.343 кг/см2, дт = 0.5 кг/см2 [5].
Рис. 7. Деформированные состояния арки (симметричная форма потери устойчивости)
На рис. 6 показана зависимость нагрузки от прогиба сродной точки арки, гдо сплошной линией обозначено решение для симметричной формы изгиба, а штриховой для несимметричной.
На рис. 7 приведено начальное и несколько деформированных состояний для симметричного изгиба, на рис. 8 для несимметричного. Численное решение выявило. что для несимметричной формы начальный этап закритического деформирования является устойчивым, и рост прогиба сопровождается ростом нагруз-
Е
Вис. 9. Зависимость деформаций от времени
ки до величины q = 0.728 кг/см2, и только затем равновесие форм становится неустойчивым. Таким образом, алгоритм, основанный на методе продолжения по параметру, показал себя в расчетах достаточно устойчивым, надежным и удобным методом решения нелинейных задач МДТТ.
6. Растяжение вязкоупругого бруса. Рассмотрен брус квадратного поперечного сечения, на одном конце которого заданы напряжения, а другой защемлен. В качестве реологической модели выбрана обобщенная модель Максвелла (11). На рис. 9 приводится сравнение точного решения (штриховая линия) и решения, полученного настоящей методикой. Как видно, результаты практически идентичны при малых деформациях и начинают расходиться лишь при увеличении величины деформации.
7. Задача изгиба конструкции, конечноэлементная модель которой представлена на рис. 10. В силу симметрии рассмотрена половина конструкции, которая разбита на две подконструкции: четырехугольное основание и обод. Примем, что левый конец основания защемлен, вертикально направленная вниз нагрузка равномерно распределена по верхней грани обода и изменяется по закону q = q0t. Задача решалась за 150 шагов при следующих параметрах: внутренний радиус обода равен 10 см, внешний 20 см, ширина левого края основания 5 см, правого 10 см, длина основания - 20 см, толщина конструкции - 1 см, E = 2.0 • 106 кг/см2, ^ = 0.3, П =1.0 • 107 кг/см2, qo = 2.0 кг/см2, aT = 1500 кг/см2. Четырехугольное основание было разбито на 20 конечных элементов по длине, 10 по ширине, 1 по толщине: обод на 20 элементов по окружному направлению, 10 по радиальному, 1 по толщине. Реологическая модель записана в виде (11), условие пластичности в виде (7). Также для наглядности была решена аналогичная задача, но без учета деформаций ползучести. На рис. 11 представлено деформированное состояние конструкции с распределением интенсивности пластических деформаций. На рис. 12 показана зависимость перемещения точки A от нагрузки, на рис. 13 -зависимость интенсивности напряжений в точке B от нагрузки, где штриховыми
Рис. 10. Копечпоэлемептпая модель «подковы»
Рис. 11. Распределение интенсивности пластических деформаций
линиями обозначены решения, полученные без учета вязкости, а сплошными с учетом последних. Из приведенных рисунков видно, что учет деформаций ползучести ведет к снижению скорости роста напряжений при возрастании деформаций.
8. Закритичоскоо поведение толстостенной трубы под воздействием моментов, приложенных с обоих концов. Физические соотношения записаны в виде (6), используется условие пластичности (7). Внешний радиус трубы 5 см, длина трубы 50 см, толщина — 1 см, Е = 2000 кг/см2, / = 0.3. При решении использовались два значения предела текучести материала: ат = 250 кг/см2 и ат = 120 кг/см2.
На рис. 14 приведена зависимость нагрузка - перемещение конца трубы, где а
ат =
250 кг/см2, штриховой — для ат = 120 кг/см2. При приделе текучести ат = 250 2
деформации. При ат = 120 кг/см2 сначала появляются пластические деформации, и только затем происходит потеря устойчивости.
На рис. 15 показано распределение интенсивности пластических деформаций для ат = 250 кг/см2, на рис. 16 - для ат = 120 кг/см2. В том и другом случае равновесное состояние закритичоского деформирования является неустойчивым. Как видно из рисунков, в трубе происходит перестройка формы и появляются вмятины, в которых наблюдаются наибольшие пластические деформации.
9. Задача деформирования грунтовой насыпи под действием собственного веса и нагружония. Грунтовый массив представляется как сплошная среда, обладающая специфическими физико-механическими свойствами. Наследовано НДС грунтовой насыпи под действием собственного веса и нагрузки, равномерно распределенной
о 100 200 300
с/, кг/см2
Гис. 12. Зависимость перемещения точки А от нагрузки
о 100 200 300
д, кг/см2
Гис. 13. Зависимость интенсивности напряжений от нагрузки в точке В
по верхней грани, нижняя грань не имеет вертикальных смещений, а боковые горизонтальных (рис. 17). Рассмотрен случай плоского деформирования. Считается. что грунт однородная среда со следующими физико-механическими свойствами: модуль деформации Е = 0.160 МПа, коэффициент бокового расширения ¡л = 0.42, сцепление С = 40 КПа, угол внутреннего трения ^ =17, плотность р = 2000 кг/м3, величина нагрузки д = 0.40 МПа. Физические соотношения записаны в виде (6), в качестве условия текучести выбран критерий Мизеса Боткина (8). Процесс деформирования был разбит на два этапа. На первом этапе находилось напряженное состояние насыпи под действием только собственного веса. На втором
4.00 8.00
М>, СМ
Рис. 14. Зависимость нагрузка перемещение конца трубы
Рис. 15. Распределение интенсивности пластических деформаций для ат = 250 кг/см2
Рис. 16. Распределение интенсивности пластических деформаций для ат = 120 кг/см2
этапе определялось НДС насыпи под действием нагрузки с учетом напряженного состояния, полученного на первом этапе.
Рис. 17. Копечпоэлемептпая модель насыпи
Рис. 18. Распределение пластических деформаций
На рис. 18 показано распределение интенсивности пластических деформаций. Из рисунка видно, что грунт в районе нагружония имеет осадку, максимальные пластические деформации возникают на склонах. Зона пластических деформаций имеет небольшую площадь, поэтому при таком нагружонии для данной формы насыпи не возникает опасных участков.
Summary
A.I. Golovanov, L.U. Sultanov. Large visco-elasto-plastic deformation of solids. Finite element, analysis of strain-stress state of visco-elasto-plastic solids with large displacements, rotations and strains is proposed. An incremental method in mixed Lagrangian-Eulerian formulation is used.
Литература
1. Голованов А.И., Султанов Л.У. Численное исследование больших упругопластиче-ских деформаций трехмерных тел // Прикладная механика. 2005. Т. 41, Л' 6. С. 36 43.
2. Повдеев A.A., Трусов П.В., Нягиии Ю.И. Большие упругопластические деформации: теория, алгоритм, приложения М.: Наука, 1986. 232 с.
3. Уилкинс M.JI. Расчет упругопластических течений // Вычислительные методы в гидродинамике. М.: Мир, 1967. С. 212 263.
4. Черных К.Ф. Нелинейная теория упругости в машиностроительных расчетах. Л.: Машиностроение, 1986. 336 с.
5. Шалалиилии В.И., Кузнецов Е.Б. Метод продолжения решения по параметру и наилучшая параметризация в прикладной математике и механике. М.: Эдиториал УРСС, 1999. 224 с.
6. МсМееклпд R.M., Rice J.R. Finite-element formulations for problems of large elastic-plastic deformation // Int. J. Solids St.uct. 1975. V. 11, No 5. P. 601 616.
7. Taylor L.M., Becker E.B. Some computational aspect of large deformation, rate-dependent. plasticity problems // Comput. Motli. Appl. Mecli. Eng. 1983. V. 41. No 3. P. 251 277.
Поступила в редакцию 20.10.05
Голованов Александр Иванович доктор физико-математических паук, профессор. проректор по научной работе и информатизации Казанского государственного университета.
E-mail: Alexandr.GolovanovQksu.ru
Султанов Ленар Усманович научный сотрудник НИИ математики и механики им. Н.Г. Чеботарева Казанского государственного университета.
E-mail: Lenar.SultanovQksu.ru
