CC BY
11
УДК 539.37
Моделирование локализации деформации в разнопрочной среде
Ольга И.Кузоватова*
*
Институт математики, Сибирский федеральный университет,
пр. Свободный 79, Красноярск, 660041, Россия
Владимир М.Садовский
Институт вычислительного моделирования СО РАН, Академгородок, 660036, Красноярск, Россия
Получена 01.03.2008, окончательный вариант 05.06.2008, принята 15.06.2008
Для исследования направлений локализации деформаций в образцах из разнопрочного материала используется специальная математическая модель, обобщающая классическую модель теории упругости. Численное решение задач проводится с помощью итерационного процесса, на каждом шаге которого решаются уравнения теории упругости с начальными напряжениями на основе метода конечных элементов.
Ключевые слова: локализация деформаций, математическая модель, теория упругости
Математическая модель
Для описания напряженно-деформированного состояния разнопрочного материала, имеющего различные пределы прочности при растяжении и сжатии, будем использовать модель сыпучей среды с пластическими связями [1]. Под действием сжимающих или растягивающих напряжений, меньших коэффициента сцепления (предела прочности связей), такая среда не деформируется. Достижение предела отвечает состоянию равновесия, в котором деформация может быть произвольной положительной величиной. Напряжения выше этого предела невозможны. Регуляризованные определяющие соотношения деформирования разнопрочной среды, учитывающие упругость частиц и связующего, приводятся к системе уравнений
Здесь а — тензор напряжений, а о — тензор сцепления, a — симметричный положительно определенный тензор коэффициентов упругости частиц, Л — параметр регуляризации, П — оператор проектирования на конус допустимых деформаций C по норме |е| = л/е : a : е. Тензор деформаций е определяется по формуле
в которой u — векторное поле перемещений, звездочка означает транспонирование.
В модели сыпучей среды с абсолютно твердыми частицами, определяющие соотношения которой получаются из (1) в пределе при a ^ ж, Л ^ 0, потенциал напряжений равен а о :
* e-mail: [email protected] t e-mail: [email protected] © Siberian Federal University. All rights reserved
(1)
2e (u) = Vu + (Vu)* ,
(2)
е+Sc (е), где Se — индикаторная функция конуса C, равная нулю на конусе и бесконечности вне его. Двойственный потенциал деформаций Sk (о — оо) выражается через индикаторную функцию сопряженного конуса K, который состоит из тензоров напряжений, образующих с тензорами из C тупые углы в смысле скалярного произведения, задаваемого сверткой о : е.
В рамках регуляризованной модели с определяющими соотношениями (1) потенциал напряжений вычисляется по формуле
ф(е) = 2(|е|2 — ТТЛ 1П(е — : ^о)12
двойственный потенциал деформаций — преобразование Юнга для Ф(е) — равен
*(о) = \( ИО + л|о — оо|0 — л|п(о — оо)|0
где п — проектор на конус К по норме |а|о = Vа : а-1 : а.
Разнопрочность материала будем описывать на основе условия Мизеса-Шлейхера, в соответствии с которым
К = {а | т(а) < жр(а)}.
Здесь т(а) = у/а' : а'/2 — интенсивность касательных напряжений, штрих означает девиа-тор тензора: а' = а + р(а) 6, р(а) = —а : 6/3 — гидростатическое давление, ж — параметр внутреннего трения, 6 — символ Кронекера.
В пространстве главных напряжений
í ^ + g2 + g3 'In П
р(а) =----, а = 0 а2 - Р 0
3
Выражение для т(а) принимает вид
Ст! - p 0 0
0 а 2 — p 0 0 0 аз — py
т(-) = </6 5> - -)2-
V »>j
Уравнение конической поверхности K записывается в форме
\J(-2 - -з)2 + (-3 - -l)2 + (-1 - -2)2 = -у/2/3ж(-1 + -2 + -з). (3)
Ось конуса совпадает с главной октаэдрической осью. Косинус угла у при вершине конуса может быть вычислен через скалярное произведение единичного вектора (-1,-2,-3), принадлежащего поверхности, и направляющего октаэдрического вектора -(1,1,1)/%/3:
cos у = —^(-1 + -2 + -з), -2 + -I + -3 = 1. (4)
3
Учитывая (3) и (4), после преобразований получим:
tg у =v/2/3®. (5)
Сопряженный конус C является круговым и его ось также совпадает с главной октаэдрической осью. Угол при вершине равен п/2 - у. Принимая уравнение сопряженного конуса в следующем общем виде
C = {е | 7(е) < ж0(e)},
где y = л/2 е' : е' — интенсивность сдвига, а 0(e) = е : S — деформация объема, определим параметр ж.
В системе главных осей тензора деформаций
/ei - 0/3 0
0(e) = ei + е2 + ез, е' = 0 е2 - 0/3
00
Выражение для y (е) принимает вид
Y (е)^/| Е(е< - )2•
V »>j
Уравнение конической поверхности конуса C приводится к форме
\/(е2 - ез)2 + (ез - е1)2 + (е1 - е2)2 = ^/3/2® (е1 + е2 + ез).
Повторяя выкладки, связанные с вычислением угла у, получим
tg (п/2 - у) = v/372 ж = ctg у,
откуда с учетом выражения (5) следует, что ж = 1/ж.
Согласно условию прочности Мизеса-Шлейхера, при достаточно больших значениях ж на некоторых элементарных площадках, проходящих через заданную точку среды, допускаются растягивающие (положительные) нормальные напряжения. Минимальное такое значение соответствует касанию конической поверхности главных координатных плоскостей отрицательного октанта. Косинус угла у, при котором происходит касание, равен скалярному произведению единичных векторов -(1,1,1)/л/3 и -(1,1, 0)/%/2, направленных по оси конуса и по одной из линий касания. Следовательно,
cos у = v/2A tg у = 1/л/2, ж = а/3/2.
В действительности растягивающие напряжения могут возникать в сыпучей среде с частицами достаточно большого размера и специфической формы, например, в кирпичной кладке. В этом случае возможно равновесие массива, находящегося в поле тяжести, с образованием нависающих свободных поверхностей. Угол внутреннего трения в такой среде превышает п/2, а в условии предельного равновесия параметр ж больше, чем л/3/2.
В случае кругового конуса Мизеса-Шлейхера соотношения (1) для вычисления напряжений по заданным деформациям в среде с упругими частицами принимают тензорный вид, и для их реализации не требуется перехода к главным осям. Для изотропной сыпучей среды тензор ао имеет вид
ао = (Ts/ж) S,
где ts — коэффициент сцепления. При вычислении проекции на конус C можно применить общую формулу (см., например, [2])
a : П (a-1 : s) = s - п (s),
которая приводит к вычислению проекции на конус Мизеса-Шлейхера по следующим формулам (k — модуль объемного сжатия, ^ — модуль сдвига):
0 0
ез - 0/3
1) если т (в) < жр (в), то а = в;
2) если т (в) > жр (в) и рр (в) + ж к т(в) < 0, то а = 0;
3) если т (в) > жр (в) и рр (в) + ж к т(в) > 0, то
жрИ г . , , ,, , рр (в) + жкт (в)
а = в + р (в) ^ *] - р (а) ^ * р (а)= р + ж2к •
Пусть П — плоская область с границей Г, состоящей из двух непересекающихся частей Ги и Гр, на первой из которых отсутствуют перемещения, а на второй задан вектор распределенной нагрузки р = (р1,р2). Задача состоит в определении векторного поля перемещений и = («1,М2) и тензорного поля напряжений, удовлетворяющих уравнениям (1), (2) для плоского деформированного состояния и уравнению равновесия в вариационной форме
JJ а : е(и) ¿П = Jр • и ¿7 для любого
п
с граничным условием
и| г =0
1 г и
для произвольного варьируемого поля перемещений и.
Эта задача сводится к двум вариационным принципам [1]. Поле перемещений минимизирует функционал
I (и> = //*<* (и)) ¿П - / р • и ¿Г
при условии и|г = 0, а поле напряжений получается в результате минимизации функционала
1 (а) = ^ Ф(а) ¿П
п
по всевозможным напряжениям, удовлетворяющим уравнению равновесия.
Для среды с абсолютно жесткими частицами вводится понятие безопасных нагрузок р, при которых среда не деформируется. Граничные точки множества безопасных нагрузок образуют множество предельных нагрузок. Для любой нагрузки р коэффициентом запаса называется число т, при котором нагрузка т • р является предельной. Справедлива следующая теорема о предельном равновесии: коэффициент запаса не превосходит кинематического коэффициента, вычисляемого по формуле:
/ ао : е(и) ¿П
т < т* = п—т---, (6)
] р • и аГ
гР
в которой и — произвольное допустимое поле перемещений, удовлетворяющее граничному условию и|г =0 и включению е(и) € С.
Локализация деформаций
В качестве примера рассмотрим плоское деформированное состояние однородного цилиндрического образца радиуса Д с радиальным разрезом, на берегах которого действует
давление р = (0,ро), ро > 0, вызванное, например, температурным расширением вставленной в надрез тонкой металлической пластины (рис. 1).
Пусть и = (иг, и2) — допустимое поле перемещений, описывающее локализацию деформации простого сдвига с дилатансией в узкой линейной зоне толщины Н, наклоненной под углом а к линии надреза (рис. 1). В декартовой системе координат, связанной с этой зоной,
иг = 7о Х2, и2 = го Х2 (0 < х2 < Н). (7)
В остальной части, вне зоны локализации, перемещения постоянны и определяются по непрерывности. После вычисления деформаций по формулам (2) из условия г (и) Е С получим
7о < ^о, ^ =\!~Ь - 4. (8)
В силу (6) оценка коэффициента запаса имеет вид
/ ип
ш* = — т£
п
г
ж / (р 1 иг + р2и2) ¿7'
гр
где ип — нормальное перемещение точки границы, точная нижняя грань берется по всевозможным перемещениям вида (7) при условии (8). Получаемая таким образом наилучшая верхняя оценка давления и наиболее вероятный угол выхода зоны локализации деформаций равны
q* = ш* • ро = — 5 =, а = arctg V. (9)
д/1 — ж2/3
Можно показать, что угол а равен углу внутреннего трения разнопрочного материала [2].
Судя по (9), в пределе при ж ^ 0 направление линии локализации перпендикулярно надрезу. При ж ^ л/3/2 зона разворачивается и становится продолжением надреза. В случае ж > л/3/2 локализации рассматриваемого типа не происходит.
Верхняя оценка (9) заведомо не является точной. Ее можно улучшить за счет введения криволинейной зоны локализации деформаций. Друккер и Прагер [3], основываясь на ассоциированном законе пластического течения, высказали гипотезу о том, что линиями локализации в плоском деформированном состоянии связной сыпучей среды могут служить только логарифмические спирали. Это следует из того, что для логарифмической спирали
угол между вектором точки и направлением касательной остается постоянным вдоль всей линии. Поэтому в узкой зоне локализации постоянным оказывается отношение нормальной и касательной составляющих вектора скорости (вектора малых перемещений), что в точности соответствует кинематике деформирования. Приведем более строгое обоснование этой гипотезы: докажем, что среди всех возможных линий локализации наименьшую верхнюю оценку дает логарифмическая спираль.
Рис. 2. Геометрические построения
Пусть Р0 — мгновенный центр вращений в плоском абсолютно жестком движении, совершаемом той частью цилиндрического образца, которая примыкает к надрезу и отделена от оставшейся неподвижной части зоной локализации деформаций (рис. 2а). Полярные координаты точки Р0 обозначим через г0 и ф0. Рассмотрим вспомогательную полярную систему координат (г, ф) с центром в точке Р0. Уравнение линии локализации можно записать в виде г = г (ф). Касательное и нормальное перемещения в слое толщины Н вблизи этой линии изменяются линейно по толщине, следовательно, деформации определяются через перемещения на верхней границе слоя по формулам
и и
£тт = О (1), епп = + О (1), 2 £„т = Н + О (1)
(10)
(п и т — нормальное и касательное направления). Имеем следующие формулы для инвариантов тензора деформаций:
в (е) = £тт + £пп, 7 (е) =
2( _ )2 + 2 2 , 2 2 + 2
3 (етт е пп ) + з е тт + з е пп + 7пт
Учитывая (10), получим
в(е) = ип + О(1), 7(-е)^3 (ЦТ)2 + (Т)2 + О(')-
Условие е(и) Е С выполняется, если ж7(е) < в(е), т. е. если
й (£)2 + (X)2 < Ж (Т + О(1)).
В пределе при Н ^ 0 получим
ит < .
а
При вращательном движении вокруг центра Po вектор малых перемещений и направлен перпендикулярно радиус-вектору точки, и его модуль вычисляется через угол поворота ф как произведение гф. Поэтому
и^ + иП = г2 ф2, ит = гф cos а, ип = гф sin а, (12)
где а — угол между вектором перемещений и касательной, зависящий, вообще говоря, от положения точки на кривой. В соответствии с представлением (12) и неравенством (11), проведем цепочку преобразований
ит • ип , ип ^ иП 1
cos а = —-, sin а = —-, tg а = — ^ - = —,
гф гф ит vun V
в итоге получим
1
tg а > -.
V
Из геометрической картины, изображенной на рис. 2б, следует, что
¿г
tg а = — .
г ау
Тогда
¿г ¿у
--= tg а ау ^ —.
г V
В результате интегрирования можно установить вид допустимых линий локализации: это только те линии локализации, на которых выполняется условие
г (ф) < гое-(^. (13)
Верхняя оценка предельного давления преобразуется к виду
Л
J мп ¿в
* Тз О 9 = —
ж
0
u y dx
Здесь Pi — точка выхода линии локализации на границу образца, ds — элемент дуги, который, в соответствии с формулами
2 2 2т2 í dr , dr
ds = dr + r d<^> , d<^> = — ctg a —, равен ds =--.
r sin a
Итак, в случае криволинейной зоны локализации
22
* Ts r2 — r2
q = —
ж Д2 +2 го Д cos у0,
Учитывая условие допустимости линии локализации (13), выберем свободные параметры в данном выражении из соображений минимума q* следующим образом:
Уо = 0, ri = roe-Vl/v,
тогда
* Ts T
q =
0 (1
. (14)
ж R2 + 2 r0 R V 7
Отсюда ясно, что логарифмическая спираль действительно обладает экстремальным свойством, хотя на самом деле точно такое же значение q* дает любая другая линия локализации с фиксированной точкой выхода на границу образца. Полярный угол этой точки может быть найден из рассмотрения рис. 2б по теореме косинусов для треугольника OP0Pi, согласно которой R2 = rQ + г2 — 2 ro г i cos р i. Расстояние до мгновенного центра вращений го также является свободным параметром в формуле для верхней оценки. Разделим числитель и знаменатель (14) на rQ и введем обозначение
Z = R,
r0
тогда
*_ Ts 1 — e-2 ^ /v
q=
ж Z2 + 2 z
Получим явное представление величины Z от угла р i:
(15)
z2
= RQ = rQ + r 2 — 2 To r i cos pi = 1+/гЛ +2 H
+ 2 — cos р i,
0 0 r0 r0
так как r i/r0 = e-^1 /v, то Z (р i) = V1 + e-2 /v + 2 e-^1 /v cos р i = R/r0. Итак, далее задача сводится к поиску минимума функции (15). Применяя асимптотическое разложение
Z2 » 1 + (1 — 2р i/v + 2pQ/v2) — 2 (1 — рi/v + pQ/(2 v))(1 — pQ/2) « (1 + 1/v2)pQ,
можно показать, что при малых р i оценка (15) совпадает с (9).
Рис. 3. Зоны локализации — логарифмические спирали
Результаты численного решения задачи минимизации для 0 < ж < л/3/2 представлены кривой 1 на рис. 3а. Для сравнения, кривая 2 на этом рисунке отвечает оценке (9). Как видно, между оценками имеется значительное расхождение, увеличивающееся с ростом ж. Дополнительный анализ показывает, что такое расхождение возникает из-за учета неоднородности поля перемещений в криволинейной зоне локализации деформаций. В частности, параметру ж = л/3/2 в обоих случаях соответствует линейная зона локализации, продолжающая разрез, однако разница в том, что оценка (9) строилась с помощью однородного по
a
длине поля перемещений, а в случае оценки (15) использовалось поле перемещений с мгновенным центром вращений в точке 0).
Серия логарифмических спиралей для различных значений параметра ж в диапазоне его изменения изображена на рис. 3б (кривые 1, 2, 3, 4 и 5 соответствуют ж = 0, 0.3, 0.5, 0.7 и 0.8).
Результаты расчетов
Численное решение задачи строится на основе метода конечных элементов [4]. Область решения разбивается на треугольные элементы. Поле перемещений ищется в кусочно-линейной форме. Вводится вектор обобщенных перемещений в вершинах треугольников. Строится глобальная матрица жесткости. В упругой задаче система линейных алгебраических уравнений решается методом сопряженных градиентов [5]. Для исследования модели разнопрочной среды применяется итерационный алгоритм метода начальных напряжений. В качестве нулевого приближения используется решение упругой задачи. Идея алгоритма состоит в замене определяющих уравнений (1) итерационной формулой (п = 1, 2, 3,...)
1
1 + Л
а : П (V
n-1
- а 1 : ао) .
Итерации продолжаются до тех пор пока норма разности двух последовательных решений не станет меньше предварительно заданной величины.
Алгоритм реализован в виде программы на алгоритмическом языке Borland C++.
Задача о деформации цилиндрического образца с надрезом. На рис. 4 а представлено поле интенсивности сдвига, полученное на основе классической теории упругости, на рис. 4 б — полученное на основе модели разнопрочной среды с параметром внутреннего трения ж = 0.3. Произведено наложение логарифмической спирали, показывающей направление локализации деформации по результатам оценочного решения.
Рис. 4. Интенсивность деформации сдвига в цилиндрическом образце
Задача о деформации прямоугольного образца с боковым надрезом. На рис. 5а
представлена схема нагружения в задаче о растяжении прямоугольного образца. Сила Р распределена по верхней и нижней границам образца равномерно. Необходимо определить векторное поле перемещений. В силу симметрии задачи была рассмотрена верхняя половина области решения. На рис. 5б представлено поле интенсивности сдвига, полученное на основе классической теории упругости, на рис. 5в - на основе модели разнопрочной среды с параметром внутреннего трения ж = 0.3. На последнем рисунке непрерывной кривой показана линия локализации деформаций - логарифмическая спираль.
а = а : £ —
а
в
Рис. 5. Схема нагружения и интенсивность деформации сдвига в образце с боковым надрезом
Задача о разрушении углеграфитового блока электролизера. Примером использования разномодульного материала может служить электролизер для получения алюминия. Электролизер является сложным высокотемпературным агрегатом (рис. 6). Основным элементом конструкции, от которого зависит срок службы, является углеграфитовый блок с жестко вмонтированным стальным блюмсом. Вследствие температурного расширения стального блюмса начинается разрушение углеграфитового блока. Обычно локализация деформаций приводит к образованию магистральных трещин. Появление подобных трещин приводит к снижению сортности получаемого алюминия и сокращению срока службы электролизера.
Рис. 6. Схема электролизера
В расчетах задавались действительные размеры поперечного сечения блоков, применяемых на Красноярском алюминиевом заводе: ширина — 0.55 м, высота — 0.4 м, ширина паза под блюмс — 0.25 м и глубина паза под блюмс — 0.15 м.
На рис. 7 а представлено поле интенсивности сдвига, полученное на основе классической теории упругости, на рис. 7б — на основе модели разнопрочной среды с параметром внутреннего трения ж = 0.3. По рис. 7б можно сделать вывод, что при рассматриваемой схеме нагружения деформации не локализуются в направлении верхней границы углеграфитово-го блока.
Рис. 7. Интенсивность деформации сдвига в блоке электролизера
Выводы
В работе на основе метода конечных элементов предложен вычислительный алгоритм для решения задач локализации деформаций в разнопрочной среде. В качестве примера рассмотрено плоское деформированное состояние однородного цилиндрического образца с радиальным разрезом и прямоугольного образца с боковым надрезом. Проведено сравнение численных результатов с точными оценочными решениями для локализованных зон деформации, которое показало, что полученные численные результаты качественно совпадают с теоретическими. В рамках модели разнопрочной среды выполнены расчеты деформации углеграфитового блока алюминиевого электролизера при тепловом расширении блюмса.
Работа выполнена при финансовой поддержке РФФИ (код проекта 08-01-00148).
Список литературы
[1] В.П.Мясников, В.М.Садовский, Вариационные принципы теории предельного равновесия разнопрочных сред, Прикладная математика и механика, 68(2004), вып. 3, 488-499.
[2] О.В.Садовская, В.М.Садовский, Математическое моделирование в задачах механики сыпучих сред, М., Физматлит, 2008.
[3] Д.Друккер, В.Прагер, Механика грунтов и пластический анализ или предельное проектирование, Определяющие законы механики грунтов, Сер. Новое в зарубежной науке, Вып. 2, М., Мир, 1975, 166-177.
[4] Г.И.Марчук, В.И.Агошков, Введение в проекционно-сеточные методы, М., Наука, 1981.
[5] Г.И.Марчук, Методы вычислительной математики, М., Наука, 1980.
Modeling of a Deformation Localization in a Medium with Different Strengths
Olga I.Kuzovatova Vladimir M.Sadovsky
A special mathematical model which is a generalization of a classical model of the elasticity theory is used for the analysis of directions of a deformation localization in samples made of solid material with different strengths. Numerical solution of the problems is carried out by means of iterative process. At each step of this process the equations of the elasticity theory with initial stresses are solved by means of the finite-elements method.
Keywords: deformations localization, mathematical model, elasticity theory
