/
Автор: Покрас И.Б.
Теги: технология обработки без снятия стружки в целом: процессы, инструмент, оборудование и приспособления математика обработка металлов металлы моделирование математическое моделирование
ISBN: 5-7526-0112-6
Год: 2002
Текст
МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ
ИЖЕВСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ
ОТДЕЛ ОБСЛУЖИВАНИЯ
УЧЕБНОЙ ЛИТЕРАТУРОЙ
И.Б. Покрас
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
ПРОЦЕССОВ ОБРАБОТКИ
МЕТАЛЛОВ ДАВЛЕНИЕМ
Рекомендовано Учебно-методическим объединением вузов по об-
разованию в области машиностроения и приборостроения в каче-
стве учебного пособия для студентов высших учебных заведений,
обучающихся по специальности 120400 «Машины и технология
обработки металлов давлением»
Ижевск 2002
УДК 621.73.01(075)
П48
Рецензенты:
А.М. Дмитриев, д-р техн, наук, пр
чл.-кор.
АНРФ, зав. каф. «Физикотехнологические процессы и оборудование
автоматизированной обработки давлением (МТ-6)»; В. Б. Дементьев,
д-р техн, наук, проф., зав. лаб. ИПН УрО РАН
И.Б. Покрас
Математическое моделирование процессов обработки ме-
П48 таллов давлением: Учеб, пособие. - Ижевск: Изд-во ИжГТУ,
2002. - 168 с.
ISBN 5-7526-0112-6
В учебном пособии изложены основные принципы математиче-
ского моделирования процессов обработки металлов давлением. Да-
на постановка соответствующих краевых задач. Приведены методы
их решении. Предназначено студентам вузов, обучающимся по спе-
циальности 120400 «Машины и технология обработки металлов дав-
лением», а также всем студентам технических вузов, специализирую-
щимся в области механики сплошной среды.
© Покрас И.Б., 2002
© Издательство ИжГТУ, 2002
НАУЧНАЯ БИБЛ И ?г*<д
ИЖЕВСКОГО ГОСУДАРСГЬьнпОГО
L таического ЖИ6ЕРСИУТД i.
ВВЕДЕНИЕ
Совершенствование процессов обработки металлов давлением
и разработка новых, высокоэффективных технологий требуют созда-
ния методов достаточно точного количественного их описания, т.е.
создания методов математического моделирования. Так как обра-
ботка металлов давлением связана с пластическим формоизменением
при определенном температурном режиме, то математическая модель
должна давать возможность определять в области пластического те-
чения распределение скоростей, деформаций, температур, напряже-
ний; рассчитывать вероятность разрушения металла и оптимизиро-
. вать условия его деформирования.
Научной основой для разработки таких моделей является теория
обработки металлов давлением, базирующаяся на основных положе-
ниях механики сплошных сред, включая теории упругости, пластич-
ности, теплопередачи, деформируемости, разрушения и др.
Для
определения напряженно-д
рмированного состояния и тем-
пературного режима приходится решать соответствующую краевую
задачу.
Наиболее универсальными для решения таких задач, с точки зре-
ния их использования для анализа пластического формоизменения
сложных геометрических объектов (поковки, штамповки, прокат
и т.д.), являются методы конечных и граничных элементов.
Основу настоящего учебного пособия составил курс лекций, чи-
таемый автором студентам специальности 1204 «Машины и техноло-
гия обработки металлов давлением» ИжГТУ. Следует отметить: по
изучаемому курсу этой специальности учебников и учебных пособий
нет, что и побудило написать эту книгу. В ней приведены способы
построения математических моделей процессов обработки металлов
давлением, позволяющие найти в области пластического течения
распределение скоростей, деформаций, температур напряжений
и рассчитать вероятность разрушения материала, а также оптималь-
ные условия его деформирования.
Даны постановка и методы решения краевых задач теории пла-
стического течения, включающие методы конечных и граничных
элементов.
3
Глава 1
ОСНОВНЫЕ ПОЛОЖЕНИЯ
МЕХАНИКИ СПЛОШНЫХ СРЕД
В этой главе, следуя [1], кратко приводятся основные формулы
теории напряжений и деформаций; при этом выделяются сведения,
наиболее важные для построения теории пластичности.
1.1. Напряженное состояние деформируемой среды
1.1.1. Напряженное состояние
В данной точке сплошной среды напряженное состояние характе-
ризуется симметричным тензором напряжения
(1.1)
где сгх, ог - нормальные, а тх/,
xxz - касательные напряжения на
площадках, перпендикулярных
к координатным осям х, z.
Вектор напряжения р на про-
извольно ориентированной пло-
щадке с единичной нормалью п
(рис. 1.1) определяется формула-
ми Коши:
Рх Wx + ^хуру
Ру ХХ)^Х ^у^у Xyz^Z^ ► О *2)
Рис. L1. Схема к выводу формул
Коши
Pz XXZ^X + Xyz”y + ^Z^Z’
где я , n , n - составляющие единичного вектора нормали и, равные
Л у л —
направляющим косинусам cos(n, х), cos (и, у)> cos (л, z).
4
Проектируя вектор р на направление нормали, получаем нор-
мальное напряжение стя, действующее на рассматриваемой площадке
стп ахпх +аупу -t-ciznz +2ххупхпу +2xyznynz у2х xznxnz.
(1.3)
Величина касательного напряжения равна
(1.4)
В каждой точке среды существуют такие три взаимно перпенди-
кулярные площадки, на которых касательные напряжения равны ну-
лю. Направления нормалей к этим площадкам образуют главные на-
правления тензора напряжения и не зависят от исходной системы ко-
ординат х, у, z. Это означает, что любое напряженное состояние
в рассматриваемой точке может быть вызвано растяжением окрест-
ности точки в трех взаимно перпендикулярных направлениях. Соот-
ветствующие напряжения называются главными нормальными напря-
жениями1, будем обозначать их через 01,02, причем условимся ну-
меровать главные оси так, что
O^Oj^Oj. (1.5) 3
Тензор напряжения, отне-
сенный к главным осям, имеет
О
О
O1 О
О с2
О о3
Нетрудно найти по форму-
лам (1.2)-(1.4), что в сечениях,
делящих пополам углы между
главными плоскостями и про-
ходящих через главные оси J,
2, 3 соответственно (рис. 1.2),
Рис. 1.2. Схема к расчету главных
касательных напряжений
касательные напряжения по абсолютной величине равны
5
Касательные напряжения в этих сечениях достигают экстремаль-
ных значений и называются главными касательными напряжениями.
Определим последние формулами:
(1-6)
С изменением ориентации площадки изменяется и величина дей-
ствующего на площадке касательного напряжения тл. Наибольшее
значение т„ в данной точке называется максимальным касательным
напряжением tmax. Если условие (1.5) выполняется, то тгпах = —т2.
Нетрудно определить по формуле (1.3), что нормальные напряже-
ния на площадках, на которых действуют главные касательные на-
пряжения (1.6), равны соответственно полусуммам
+ с3
°3 +<т1
°1 +о2
(1.7)
2
9
2
Тензор напряжения можно задать, указав главные напряжения Oi,
аг, стз и главные направления 1, 2, 3. Такое задание отличается меха-
нической ясностью. Главные напряжения ст,, (i = 1, 2, 3) являются кор-
нями кубического уравнения
или
-А.3+7,(ТО)Х2 +/2(ТО)Х+ 73(Та) = 0.
(1-8)
Нормальное напряжение стл на заданной площадке не зависит,
очевидно, от выбора координатной системы и изменяется лишь при
повороте площадки. Главные напряжения сть Стг, Оз являются экстре-
мальными значениями нормального напряжения ст„ и также не зави-
сят от выбора координатной системы. Уравнение (1.8) может быть
получено как условие экстремума стп. Следовательно, коэффициенты
кубического уравнения (1.8) не изменяются при переходе от одной
прямоугольной системы координат к другой, т.е. инвариантны. Эти
коэффициенты
6
/1(Та) = О1 + СТ2+СТ3 =ЗСГ,
Л(ТСТ) — + °2а3 + CT3CTl)i ”
/3(ТСТ) = 0102^3,
(1.9)
записанные, для кратности, в главных осях, называются, соответст-
венно, линейным, квадратичным и кубическим инвариантами тензо-
ра; ими удобно оперировать, так как они являются целыми рацио-
нальными и притом симметричными (т.е. не изменяющимися при пе-
рестановке аргументов) функциями компонент напряжения.
Величина ст = —(ст
называется средним (или гидроста-
тическим) давлением в точке. Смысл остальных инвариантов выяс-
нится ниже.
1.1.2. Девиатор напряжения
Так как материалы обладают, как правило, различными механи-
ческими свойствами по отношению к сдвигу и равномерному всесто-
роннему сжатию, выгодно представить тензор напряжения в виде
суммы:
T^oTj+Д,, (1.10)
где стТ, = О
О О
ст О
О ст
- шаровой тензор, соответствующий среднему
О
давлению в точке;
(1-11)
- тензор, характеризующий касательные напряжения в данной точке
и называемый девиатором напряжения; Tj - так называемый единич-
ный тензор
7
для которого любое направление является главным, а диагональные
элементы в произвольной прямоугольной системе координат х, у, z
равны единицам.
Нормальные составляющие последнего (т.е. о, а - с, ст,- а) бу-
дем иногда обозначать через S S_. Главные направления девиа-
тора напряжения D0 и тензора напряжения То совпадают, а главные
значения S{ отличаются от а, на величину среднего давления и опре-
деляются, очевидно, кубическим уравнением
-K+I2(D^ + I3(Da)=Q,
(1.12)
все корни которого также вещественны.
Инварианты девиатора легко получить из (1.9), если заменить о,,
о2, ст3 соответственно на S2, S3 :
pi(Z>e) = O;
Очевидно, что девиатор напряжения характеризуется лишь пятью
независимыми величинами.
Неотрицательную величину
(1.14)
называют интенсивностью касательных напряжений.
В ряде случаев используют также величину сти, называемую интен-
сивностью напряжений
ои = Тэт.
Интенсивность касательных напряжений обращается в нуль толь-
ко в том случае, когда напряженное состояние является состоянием
гидростатического давления.
Для чистого сдвига aj = т, <з2 = 0, о3 = -т, где т - напряжение сдвига.
Следовательно, Т = т.
В случае простого растяжения (сжатия) в направлении оси х
= *i; Gy = °г = тхУ = тогда
8
(1.15)
Поскольку кубическое уравнение (1.12) имеет вещественные кор-
ни, решение его находится в тригонометрическом виде. Используя
известные формулы алгебры, можно выразить главные компоненты
девиатора через инварианты
(1.16)
Угол со», определяется из уравнения
. зЛг3(па)
- cos3o>- =---
2Т3
(1.17)
Нетрудно найти, исходя из (1.6), (1.16), главные касательные на-
пряжения:
(1.18)
Угол сПз, как известно, изменяется в пределах
(1-19)
Действительно, т.к. aj > с2 ст3, то -ц > 0, т2 > 0, т3 > 0, т.е.
sin(oa > 0, откуда следует (1.19).
Как уже указывалось, тта, = -т2- Отсюда вследствие (1.19) вытека-
ет неравенство, установленное иным путем А.А. Ильюшиным:
(1.20)
9
Таким образом, интенсивность касательных напряжений Т и мак-
симальное касательное напряжение тгоат незначительно отличаются
друг от друга. Так,
7^1,08^ (1.21)
с наибольшей погрешностью около 7 %.
1.1.3. Тензорные обозначения
При рассмотрении общих вопросов теории пластичности исполь-
зование тензорных обозначений упрощает изложение и делает его
более ясным. Тензорные обозначения находят все более широкое
распространение в современной научной литературе по теории пла-
стичности. По отмеченным причинам тензорные обозначения ис-
пользуются и в ряде разделов этой книги.
Декартовы координаты х, у, z будем обозначать через xt, х2, хз
и записывать их как хр где индекс i принимает значения 1, 2, 3. Разу-
меется, вместо i можно взять другую букву (например J, j ~ 1, 2, 3;
обычно используются латинские буквы). Через п (или, скажем, nJ)
обозначим составляющие единичного вектора нормали к площадке;
очевидно, что равны направляющим косинусам нормали.
Компоненты тензора напряжения можно теперь обозначить через
i, j = 1, 2, 3. Вследствие закона взаимности касательных напряжений
а.. = о... Соотношения между тензорными обозначениями и использо-
ванными выше «техническими» обозначениями очевидны: о^ст*,
012 -хху и т.д. Условимся далее говорить о тензоре напряжения как
о тензоре с,,.
Формулы Коши (1.2) можно теперь представить в форме
з
Р t = J~ 1» 2, 3.
i=l
Широкое распространение получило правило суммирования, вве-
денное А. Эйнштейном. Опустим знак суммы, принимая, что по вся-
кому дважды повторяющемуся в одночлене латинскому индексу про-
водится суммирование по значениям 1, 2, 3. Тогда предыдущие фор-
мулы переписываются в виде
Pj=^ijni- (1-22)
10
Повторяющийся индекс i называется немым (или индексом сумми-
рования); он может быть заменен любой другой (обычно латинской)
буквой. В каждом одночлене один и тот же немой индекс не должен
встречаться более двух раз. Индекс j иногда называется свободным.
Легко видеть, что нормальное напряжение стя (1.3) равно
(1.23)
здесь два немых индекса i, j, и, следовательно, проводится двойное
суммирование; свободных индексов нет.
Среднее давление равно
8,у =
Символ Кронекера (дельта-символ) определяется соотношениями
1 при z = 7,
О при i * j.
Тензор с такими составляющими в системе координат х, называ-
ется единичным тензором (см. Ti).
Девиатор напряжения имеет составляющие
Линейный вариант девиатора равен нулю, т.е. Sa = 0. Нетрудно
убедиться в том, что интенсивность касательных напряжений в но-
вых обозначениях равна
2
Среднее давление можно также записать в форме о =—ст
1.2. Деформация
1.2.1. Тензор деформации
Пусть при деформации среды точки последней получили смеще-
ние и, составляющие которого обозначим через их, и , и,. Деформация
среды характеризуется симметричным тензором деформации:
11
е
составляющие которого равны
1
f xz
2
1
Y yz
2
Г ди' дих диу ди Qu ди2
xy dy dx dx dy dx dy dx dy
(1.26)
Тензор деформации, как и всякий симметричный тензор, приво-
дится к главным осям:
£] О О
О е2 О
О О Е3
причем £], е2, е3 называются главными удлинениями. Это означает, что
всякая деформация может быть осуществлена простыми растяжения-
ми в трех взаимно перпендикулярных направлениях (главных на-
правлениях).
Разности
Y1 =е2 -83,у2 =е3-£1эу3 =Ej-е2 (1-27)
называются главными сдвигами. Наибольший по величине сдвиг
в данной точке будем называть максимальным сдвигом у тм .
1.2.2. Малая деформация
В случае малой деформации компоненты ех, е?, ..., yxz малы по
сравнению с единицей; если, кроме того, достаточно малы углы по-
ворота, то в формулах (1.26) можно пренебрегать произведениями
' дих )2 дих дих
—— , ———..., следовательно,
k дх J дх ду
12
УХ2
(128)
Если ех, е , е/ представляют собой относительные удлинения соот-
ветственно в направлениях осей х, у, z, а у , у , уХ1 - относительные
сдвиги (у - изменение угла между осями х, у, z и т.д.), относительное
изменение объема равно
е = Ех+£>,+Ег. (1.29)
Эти простые формулы непригодны, если необходимо описать зна-
чительные формоизменения массивных тел; тогда компоненты де-
формации сравнимы по величине с единицей и нужно исходить из
общих зависимостей (1.26). Подчеркнем также, что даже при малых
удлинениях и сдвигах линейные соотношения (1.28) часто оказыва-
ются недостаточными в вопросах деформации и устойчивости гиб-
ких тел (стержни, пластины, оболочки) вследствие того, что элемен-
ты тела испытывают значительные перемещения и повороты.
В дальнейшем, говоря о малой деформации, мы будем подразумевать
такую деформацию, когда формулы (1.28) применимы.
Ниже нередко используются тензорные обозначения компонент
деформации
(1.30)
где и{- составляющие вектора перемещения; xt - декартовы коорди-
наты. Легко видеть, что е = е1?6„ .
1.2.3. Инварианты
Инварианты тензора деформации образуются так же, как для тен-
зора напряжения, и в главных осях имеют вид:
/1(Т£)=е1+е2+е3; 1
^2 1Те )=—(fiiEj +Е2Ез +E3Et );> (1’31)
Л(Те)=€]Е2Е3. J
Удобно пользоваться представлением тензора деформации в виде
суммы
T^I/ЗеТх + Л, (1.32)
13
где 1/ЗеТ] - шаровой тензор, соответствующий объемному расшире-
нию, а девиатор деформации
характеризует изменение формы элемента среды, обусловливаемое
сдвигами. Инварианты девиатора деформации равны
Л(Ре) = 0;
(1.33)
В теории пластичности важную роль играет квадратичный инва-
риант /2(£>£), который можно рассматривать как суммарную характе-
ристику искажения формы элемента среды. Неотрицательная вели-
чина называется интенсивностью деформаций сдвига.
В механике сплошных сред используют и другую величину, назы-
ваемую интенсивностью деформаций
Численный множитель перед корнем в (1.34) выбран так, чтобы
при чистом сдвиге интенсивность Г равнялась величине сдвига у.
Соотношение (1.32) может быть записано также в форме
(1.35)
14
где etj - компоненты девиатора деформации. Первое равенство (1.33)
в этих обозначениях имеет вид etj = 0, а интенсивность деформаций
сдвига равна
Г = (2е,7еД (1.36)
1.2.4. Условия совместности деформаций
Компоненты деформации должны удовлетворять шести тождест-
венным соотношениям Сен-Венана:
а2ех Э2е> a2YXj,
ду2 дх2 дхду
& xz ху
дудг дх дх ду dz )
(1.37)
Остальные соотношения получаются из выписанных круговой за-
меной индексов.
1.2.5. Компоненты деформации в цилиндрических
и сферических координатах
В дальнейшем нам понадобятся выражения компонент деформа-
ции в цилиндрических и сферических координатах; приводим их без
вывода.
Цилиндрические координаты г, ф, г. Пусть компоненты вектора
смещения иг, и^, иг не зависят от ф; тогда относительные удлинения
и сдвиги имеют вид:
(1-38)
Сферические координаты г, ф, £. В интересующем нас случае цен-
тральной симметрии компоненты вектора смещения uv = мх = 0, а
~ &Г’ еФ=ех=~» = Тфх ~ 7а = 0- (1-39)
15
1.3. Скорость деформации
1.3.1. Тензор скорости деформации
Пусть частицы среды движутся со скоростью о, составляющие ко-
торой равны их = иДх, у, z, t), vy = иу(х, у, z, t), и, = иг(х, у, z, t).
В течение бесконечно малого промежутка времени dt среда испы-
тывает бесконечно малую деформацию, определяемую перемещениями
vjdt, v2dt. Компоненты этой деформации, вычисленные по (1.28),
имеют общий множитель dt, разделив на который, получаем компо-
ненты симметричного тензора скорости деформации
где
Величины Е , определяют скорости относительных удлинений
элементарного объема в направлениях координатных осей; т]х?, т|>2, г|х2
определяют угловые скорости скашивания первоначально прямых уг-
лов. Скорость относительного объемного расширения равна
^ = ^+^+^-=divu-
(1.41)
Кроме скорости чистой деформации, характеризуемой тензором
Д, элементарный объем испытывает жесткое смещение, определяе-
мое поступательной скоростью и и вращением с угловой скоро-
стью со = l/2rot и.
Ускорение движущейся частицы среды определяется полной
(субстанциональной) производной скорости
(1 -42)
16
Здесь первый член характеризует локальные изменения, ос-
тальные - представляют трансляционную часть, учитывающую
изменения вследствие переноса частицы в соседнюю точку про-
странства.
В тензорных обозначениях компоненты скорости деформации
равны
где о. - компоненты вектора скорости.
1.3.2. Инварианты тензора скорости деформации
Инварианты тензора и девиатора можно получить из фор-
мул (1.31), (1.33) заменой ех,..., уЛ2 на 2,Л,..., *]„. Выпишем лишь выра-
жение интенсивности скоростей деформации сдвига:
_________________________
= + & -^х)2 + л2лг + nt-)- (1.43)
Введем также в рассмотрение величину £и, называемую интенсив-
ностью скоростей деформаций:
Сумма интенсивностей последовательных малых деформаций, ко-
торую претерпела рассматриваемая частица с момента возникнове-
ния в ней пластических деформаций (т = 0) до данного момента т — t,
называется степенью деформации и определяется как
t
(1.45)
о
Аналогичным образом определяется степень деформации сдвига:
t
Л=[ЯА. (1.46)
о
Интегралы в (1.45) и (1.46) должны быть подсчитаны, естественно,
вдоль траекторий движения частиц.
НАУЧНАЯ БИБЛ ИОТЕКА
ИЖЕВСКОГО ГОСУДАРСТВЕННОГО
•' ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА (
17
Пусть х, у, z - координаты движущейся по линии тока частицы,
тогда
(1*47)
dx dy dz .
— - — = — = dt.
Dj,
Зная в установившемся движении их, и , иг как некоторые функ-
ции координат, можно решить эту систему дифференцированных урав-
нений и определить линию тока любой частицы деформируемого тела.
В ряде работ компоненты тензора скорости деформации обозна-
чаются так же, как компоненты тензора деформации, но с точкой на-
верху:
где
Скорость относительного изменения объема
ёо =ёх + + ez =divu.
В случае несжимаемой среды £о = 0.
В тензорных обозначениях
(1.49)
В настоящей работе используются и те и другие обозначения.
18
Интенсивность скоростей деформаций
(1.50)
1.3.3. Деформация и скорость деформации
Так как скорости суть полные производные по времени от пере-
dU:
мещений и,- =—‘-.то
dt
if Г du Г
^=7 • О-51)
21 oXj dt ex; dt '
Очевидно, что
5.7’4^ (1-52)
В случае малой деформации существуют простые соотношения
между компонентами деформации и компонентами скорости дефор-
д
мации, именно так, как теперь и,- =—м.-, то
dt
0-53)
J dt J
\Г Ж ^ui
Ускорения определяются формулами со. =—
dt
Трансляционная часть в выражении полной производной отбра-
сывается на том основании, что при малой деформации производные
по координатам от смешения и скоростей обычно можно считать ма-
лыми.
Л. к д
Следует, наконец, отметить, что , т.к. главные оси тензо-
dt
ра деформации и скоростей деформации, вообще говоря, не совпа-
дают.
1.3.4. Приращения компонент деформации
Механические свойства металлов в условиях сравнительно мед-
ленной пластической деформации при невысокой температуре прак-
тически не зависят, как будет выяснено ниже, от скорости деформиро-
вания. В этом случае представляют интерес, по сути дела, не скорости
19
деформации, а бесконечно малые приращения dt (условно обозна-
чаем их через tfe.., помня, что эти величины не являются, вообще го-
воря, дифференциалами компонент деформации), которые опреде-
ляются,согласно (1.5 Ц формулой
(1.54)
образуют тензор и имеют простой физический смысл. Соотно-
шения (1.54) применимы для описания больших де
рмаций, кото-
рые можно получить суммированием значений бесконечно малых
изменений (1.54).
В формулах (1.54) приращения компонент деформации вычисля-
ются по отношению к текущему (мгновенному) состоянию; система
координат х( предполагается «вмороженной» в данный элементарный
объем.
Рассмотрим, например, однородное растяжение цилиндра вдоль
. dl ,
его оси, совпадающей с осью xt; тогда ае1 =—, где I - текущая длина
цилиндра; dl - бесконечно малое ее изменение.
Суммирование приводит к так называемому натуральному удли-
fdl I
нению I — = 1п-—, где Zo - начальная длина.
/о i !о
Если главные оси при деформации не поворачиваются, интегралы
J de, имеют простой физический смысл, равняясь соответствующим на-
туральным удлинениям In Z .• /10. Очевидно, что при этом справедлив
простой закон сложения деформаций: сумма последовательных нату-
ральных удлинений равна суммарному натуральному удлинению.
В общем случае интегралы ] de^ не вычисляются и не имеют оп-
ределенного физического смысла, т.е. известны компоненты de*.
в функции некоторого параметра (например, нагрузки). Это ограни-
чивает область применения натуральных удлинений как меры де-
формации случаем фиксированных главных направлений.
Инварианты тензора ТА (девиатора £)Л), которые получаются
из соответствующих инвариантов тензора Те переходом к компонен-
там de*, будем обозначать через de, dT, .
20
Подчеркнем еще раз, что величины не следует рассматривать
как дифференциалы компонент деформации £„. Это будет верно
только для малой деформации, когда справедливы формулы (1.28),
тогда имеет место простая суперпозиция деформаций, и интегралы
Г de/j суть компоненты деформации.
1.3.5. Условия совместности скоростей деформации
Компоненты скорости деформации, как и компоненты деформа-
ции, не могут быть заданы произвольно. Они должны удовлетворять
шести условиям совместности, вполне аналогичным условиям совме-
стности Сен-Венана (1.37):
1.3.6. Случай несжимаемой среды
Для несжимаемой среды Е, = 0, т.е.
^ = 0.
дх/
(1-56)
При этом условии компоненты являются компонентами девиа-
тора скорости деформации и интенсивность скоростей деформаций
сдвига равна
н = (21;,Г .
(1.57)
1.4. Дифференциальные уравнения равновесия.
Граничные и начальные условия
Дифференциальные уравнения равновесия в прямоугольной сис-
теме координат
21
dx ду dz
дх ду dz
^xz , ! S<*z
дх dy dz
(1.58)
В тензорных обозначениях эти уравнения записываются в форме
tog
дх-.
= 0,
(1.59)
или = 0.
Подстрочный индекс после запятой означает частную производ-
ную по величине, которую представляет подстрочный индекс.
В дальнейшем нам понадобятся дифференциальные уравнения
равновесия в цилиндрических и сферических координатах; приводим
эти уравнения без вывода.
Уравнения равновесия в цилиндрических координатах. В цилиндри-
ческих координатах г, <р, г уравнения равновесия имеют вид:
Эст, dz „ Л
—+------—+——+---------=0;
dr г д<р dz г
J . 2Ч +nF _Л.
---+------+----+----- + р/'ф -V,
dr г Sep dz г
дтп 1 даг
+----+ + +pF2 =0.
dr---------------г dtp dz-r
(1.60)
Уравнения равновесия в сферических координатах. В сферических
координатах г (радиус), <р (долгота), х (широта) в случае центральной
симметрии уравнение равновесия имеет вид
Г ( ) Z7 Л
аг г
(1.61)
причем стф=оя, = 0.
Граничные условия. Кроме приведенных выше уравнений, мы рас-
полагаем еще граничными условиями, которые могут иметь разнооб-
разный характер.
22
На границе S тела могут быть заданы нагрузки рх, ру, рг. В этом
случае на S должны выполняться уравнения, которые будут условия-
ми равновесия элементарного тетраэдра, примыкающего к границе,
под действием внутренних и внешних сил.
Могут быть заданы смещения (или скорости) точек границы тела.
Наконец, встречаются смешанные граничные условия, когда на
границе частично заданы нагрузки, частично - смещения (или ско-
рости).
Начальные условия. Если процесс деформации является нестацио-
нарным и описывается уравнениями, содержащими производные по
времени (или по параметру нагрузки), необходимо задать состояние
тела в начальный момент времени (начальные условия).
1.5. Условия текучести
Условия текучести представляют собой критерии перехода из уп-
ругого состояния в пластическое. Эти критерии основаны на допу-
щении, что состояние текучести может быть выражено через главные
нормальные напряжения или через компоненты напряжения и не за-
висят от пути нагружения в упругом состоянии. Поэтому условия те-
кучести могут быть выражены соотношениями
/(ст,) = О или/(ст;у)=0.
(1.62)
Функции f, входящие в состав соотношения, сохраняют свой вид
при любых напряженных состояниях.
Этим условиям текучести для рассматриваемой точки соответст-
вуют некоторые поверхности в трехмерном пространстве главных
нормальных напряжений или некоторые гиперповерхности в шес-
тимерном пространстве компонент напряжений ст.
Для простейших напряженных состояний простого растяжения
и простого сдвига условие текучести для состояния идеальной пла-
стичности может быть представлено в виде
а = а5 и t = tv.
(1.63)
Для объемного напряженного состояния наибольшее распростра-
нение получило условие текучести Мизеса-Губера, согласно кото-
рому при пластическом состоянии материала интенсивность напря-
жений сдвига (интенсивность напряжений) есть величина посто-
янная.
23
В главных напряжениях это условие имеет вид
(ci -<т2)2 +(о2 -а3)2 + (а3 -а,)2 = 2а3 . (1.64)
В компонентах тензора напряжений для прямоугольной системы
координат оно принимает вид
+ (°z 6 V ху ^ух ^xz ) — • (1.65)
Условие Мизеса-Губера может быть записано в форме
Т = или сти = и*. (1.66)
При наличии упрочнения условие текучести принимает вид
Т = ф(Г) или Т = <р(Я), (1.67)
где <р — некоторая функция упрочнения.
1.6. Соотношения между напряжениями,
деформациями и скоростями деформации
1.6.1. Простейшие реологические модели
Результаты экспериментов по растяжению цилиндрических об-
разцов позволяют установить следующие фундаментальные свойства
реального металла: упругость, вязкость и пластичность. Особенности
поведения сплошной среды под действием приложенной нагрузки
могут быть иллюстрированы комбинацией этих фундаментальных
свойств.
В связи с этим удобно ввести простые реологические модели, опи-
сывающие поведение некоторых идеализированных сред, условно
изображая их механическими элементами.
Следуя этому, рассмотрим линейное напряженное состояние
(растяжения стержня). Обозначим: а - соответствующее напряжение,
£ - относительное удлинение, £ = dddt - скорость относительного уд-
линения.
Модель линейно-упругой среды, подчиняющейся закону Гука
g = Ee, (1.68)
будем условно изображать в виде пружины (рис. 1.3).
24
Модель линейно-вязкой среды, следующей закону вязкости Ньютона:
, </е
И ~г.
(1-69)
может быть представлена в виде поршня, перемещающегося в цилин-
дре, наполненном вязкой жидкостью. При этом жидкость вытекает
через зазор между стенкой цилиндра и поршнем (рис. 1.4).
При построении модели жесткопластической среды будем пред-
полагать, что при напряжениях ниже предела текучести деформации
отсутствуют. Пластическое течение имеет место при напряжении,
удовлетворяющем условию текучести
с = (1.70)
Представим эту модель в виде груза, покоящегося на плоскости
(элемент сухого трения, рис. 1.5, а).
Перейдем к рассмотрению простейших комбинированных моде-
лей. Соединим последовательно упругий и пластический элементы
(рис. 1.5, б). В результате получим модель упругопластической среды.
Диаграмма о - е для этой среды показана на рис. 1.5, б. Общая де-
формация при этом состоит из двух частей: упругой ее и пластиче-
ской е?:
г^еР + еР. (1.71)
При снятии нагрузки упругая деформация исчезает, остается пла-
стическая деформация. На рис. 1.5, в, г представлены диаграммы о - е
для жесткопластической и упругопластической линейно упрочняю-
щейся среды.
Рис. 1.4. Модель линейно-вязкой
среды
Рис. 1.3. Реологическая модель
линейно-упругой среды
25
Рис. 1.5. Модели пластических сред: а - жесткоидеальнопластическая среда; б - уп-
ругоидеальнопластическая среда; в - жесткопластическая линейноупрочняющаяся среда;
г - упругопластическая линейноупрочняющаяся среда
Соединим последовательно упругий и вязкий элементы (рис. 1.6).
Скорость деформации = dddt есть сумма упругой составляющей
£е=—и вязкой составляющей ел’=о/ц', отвечающей одному
Е dt
и тому же напряжению:
_ 1 о
(1.72)
Это уравнение соответствует модели упруговязкой среды Мак-
свелла. Рассмотрим некоторые свойства этой среды. Пусть напряже-
ния постоянны (о = const). Тогда dcrfdt = 0 и материал течет подобно
вязкой жидкости.
Зафиксируем теперь деформацию, приложив в момент t = 0 напряже-
ние а (0) и закрепив концы стержня. Поскольку при этом dddt = 0, урав-
нение (1.72) принимает вид — —— + — = 0, откуда
Е dt р.
а = ст(О)ехр[-Г/Го], (1.73)
где величина Го = называется временем релаксации. Она пред-
ставляет собой время, за которое начальное напряжение уменьша-
ется в е = 2,718 раз.
26
Таким образом, модель среды Максвелла позволяет описать важ-
ное свойство реальных тел, заключающееся в падении по экспонен-
циальному закону напряжений при неизменной деформации (так на-
зываемая релаксация напряжении).
Перейдем к анализу параллельного соединения упругого и вяз-
кого элементов (модель упруговязкой среды Фойгта) (рис. 1.7).
Очевидно, напряжение о есть сумма упругой ае = Ее и вязкой
ctv =|i' de)dt составляющих:
о = Ее + р.
, de
~dt
(1-74)
В отличие от выражения (1.72) уравнение (1.74) не описывает про-
цесса релаксации, так как при е = const напряжение также • остается
постоянным, а среда ведет себя как упругая. Если же напряжение а
остается постоянным, то деформация постепенно нарастает по закону
(1-75)
стремясь к значению <з!Е, т.е. возникает ползучесть.
Рассмотрим комбинацию вязких и пластических свойств. После-
довательное соединение двух элементов - вязкого и пластического -
приводит к модели вязкопластической среды. Она обладает свойст-
вами линейно-вязкой среды при ст < ст’ и течет подобно идеально
пластическому телу при ст = ст,.
ц'
Рис. 1.6. Упругая среда Максвелла Рис. 1.7. Упруговязкая среда Фойгта
Параллельное соединение вязкого и пластического элементов
также дает вязкопластическую среду {среда Шведова-Бингама). Пове-
дение среды описывается уравнением
ст = ст, + p’dedt при ст > CTf;
(1-76)
при ст < стЛ деформация отсутствует.
27
1.6.2. Соотношения между напряжениями и деформациями
для упругого состояния металлов
Согласно закону Гука для упругого состояния изотропного ме-
талла, компоненты девиаторов деформаций и напряжений пропор-
циональны
1 , \ z х
ех е — с) (х, у, z),
(МЪ
где (х, у, z) означает, что остальные четыре соотношения получаются
круговой перестановкой индексов; G - модуль сдвига, связанный
с модулем Юнга Е и коэффициентом Пуассона соотношением
Е = 2G(1 + v);
Модуль объемного сжатия
(1-78)
В тензорных обозначениях уравнения (1.77) имеют вид
Е-. = —о*.
У 2G V’
где Е'у и Су - компоненты девиаторов деформаций и напряжений со-
ответственно.
Компоненты тензоров деформаций и напряжений связаны соот-
ношением
(1.79)
где б,у = 1 при i ~j; 5. = 0 при i *J.
Разрешив эти соотношения относительно компонент напряжения,
получим иную форму закона:
ст1?. = Хе8у +2це,7,
где Л и ц = G - упругие константы Ламе.
28
1.6.3. Соотношения между напряжениями и малыми упру-
гопластическими деформациями (деформационная теория)
При выводе основных соотношений деформационной теории пла-
стичности приняты следующие гипотезы:
1) тело изотропно;
2) относительное изменение объема является упругой деформацией,
пропорциональной среднему давлению:
£=кс, (1.80)
. 1 —2v
где к =----, v - коэффициент Пуассона;
Е
3) девиаторы деформации и напряжения пропорциональны
PE = W>O. (1.81)
В прямоугольной системе координат соотношения между компо-
нентами девиаторов и напряжений примут вид
ел - е = у(аЛ ~ Уху = Ъ^ху 5'
еу - е = - a); y>z = 2^yz\ ►
Ег - Е = V(<T* - о); Yzx = 2^гх-
(1.82)
Если в этих уравнениях заменить е его значением из (1.80), то по-
лучим систему уравнений пластичности для малых деформаций при
простом нагружении, известную под именем уравнений Генки:
= VI ^х
еу - <|/ ау
^2 = V| ° 2
(1.83)
У ХУ ~ 2V^xy»
У yz k
Kx = 2утгх.
Из уравнений (1.82) получим зависимость, связывающую ТиГ
(еи и <ти)
Г = 2уГ;
(1-84)
В тензорных обозначениях соотношение (1.84) можно записать
в виде
Еу = 1/ ЗАхэ8(у + • 0 -85)
29
Соотношение (1.85) нетрудно также разрешить и относительно
напряжений
(186)
к у
Неизвестную функцию у выразим через Г и Г, а также стн и е*.
В результате получим
3 еи
2 °и
(1.87)
Подставляя найденные значения у в (1.82), получим
(1.88)
или
(1.89)
Во многих случаях при решении задач теории пластичности мож-
но пренебречь сжимаемостью материала, т.е. положить v = 0,5, тогда
к — 0 и, согласно (1.80), е = 0.
Для материала, удовлетворяющего условию пластичности в фор-
ме Мизеса (Т = т,), система уравнений (1.89) примет вид
(1.90)
В случае условия текучести с упрочнением необходимо задать
значение Т - ЦГ).
30
1.6.4. Соотношения между напряжениями и скоростями
деформации (теория течения)
Общие соотношения. Процесс пластической деформации является
необратимым, большая часть работы деформации переходит в тепло.
Напряжения в конечном состоянии зависят от пути деформирования.
В связи с этим уравнения, описывающие пластическую деформацию,
в принципе не могут быть конечными соотношениями, связывающи-
ми компоненты напряжения и деформации (аналогично соотношени-
ям закона Гука), а должны быть дифференциальными (и притом не-
интегрируемыми) зависимостями.
Уравнения теории пластического течения устанавливают связь
между бесконечно малыми приращениями деформаций и напряже-
ний, самими напряжениями и некоторыми параметрами пластиче-
ского состояния.
Рассмотрим исходные положения этой теории:
1) тело изотропно;
2) относительное изменение объема мало и является упругой де-
формацией, пропорциональной среднему давлению:
е = А:о (1.91)
или tfe = kdc;
3) полные приращения составляющих деформаций cfe,-, склады-
ваются из приращений составляющих упругой деформации dzg и
пластической деформации defi
dEy = de- + defj.
(1.92)
Приращения составляющих упругой деформации связаны с при-
ращениями составляющих напряжения законом Гука
(1.93)
4) девиатор напряжения Da и девиатор приращений пластической
деформации D$e пропорциональны, т.е.
(1.94)
где dk - некоторый бесконечно малый скалярный множитель. Это
положение обобщает результаты опытов по сложному нагружению,
в которых направления главных осей и соотношения главных напря-
31
жений изменялись. Согласно экспериментам приращения состав-
ляющих пластической деформации («скорости пластической дефор-
мации») пропорциональны напряжениям в данный момент времени.
Другими словами, напряженное состояние определяет мгновенные
приращения компонент пластической деформации.
Из (1.94) вытекают соотношения (т.к. dep = 0):
— dkSjj.
(1-95)
Согласно (1-92) получаем полные приращения компонент дефор-
мации:
flfejу =<fe?+ dkSfj,
(1.96)
где приращения компонент упругой деформации следует взять со-
гласно закону Гука (1.93). В результате получим
(1.97)
Эти уравнения носят название уравнений Прандтля-Рейсса.
Если в уравнениях (1.97) пренебречь упругой деформацией, т.е.
принять схему жесткопластического тела и разделить правые и левые
части на dt, то получим систему уравнений пластичности Леви-Мизеса
е.*. — А.,..,
(1.98)
или в компонентах тензоров скоростей деформаций и напряжений
(1-99)
При условии пластичности Мизеса коэффициент X определяется
следующей формулой
(1.100)
Уравнения Леви-Мизеса относительно напряжений в этом случае
примут вид:
32
где Н - интенсивность скоростей деформации сдвига.
3 &
В общем случае Л =---— и уравнения (1.99) примут вид:
2
(1.102)
(1.101)
Рассмотрим некоторое обобщение теории пластических течений,
в котором условие текучести усложнено, а упругие деформации опу-
щены [1].
Введем в рассмотрение функцию flaX
Если эта функция обладает следующим свойством
(1.103)
то ее называют пластическим потенциалом, а выражение (1.103) - ас-
социированным законом пластического течения, поскольку последнее свя-
зывается (ассоциируется) с условием текучести. В этом соотношении
= (1.104)
2 аи
В этом случае легко устанавливаются теорема единственности
и экстремальные принципы, что сообщает теории законченность.
1.6.5. Соотношение между напряжениями и скоростями
Для вывода соотношений между напряжениями и скоростями де-
формаций примем некоторые допущения. Прежде всего будем счи-
тать, что среда несжимаема:
33
= 0.
Далее компоненты напряжения складываются из компонент
напряжений а'у, связанных с пластическими свойствами, и компо-
нент напряжений о-, вызываемых вязким сопротивлением:
° 7 &ij + °t/-
(1.105)
Компоненты напряжений с'.-определяются уравнениями теории
пластичности Мизеса, т.е.
(ЫОб)
н
причем удовлетворяется условие текучести Мизеса
5^= 2г,. (1.107)
Вязкое сопротивление определяется законом вязкости Ньютона
— 2р£у,
(1.108)
где ц - коэффициент динамической вязкости.
Складывая напряжения, приходим к соотношениям вязко-
пластической среды
(1.109)
Отсюда вытекает зависимость
Т’ = т5+рЯ.
(1.110)
1.7. Некоторые обобщения
Используя гипотезу «единой кривой», можно обобщить получен-
ные результаты и представить соотношения между напряжениями
и деформациями (или скоростями деформаций) в виде
Е,у = \|/5; или £у = .
(1.111)
34
Если у = const = 1/2G, мы приходим к соотношениям для упругой
среды (к закону Гука).
При у = Г/2т, получаем соотношения для идеально пластической
среды с условием текучести Мизеса;
при ц» = l/2g(r) ц/ - 1/2 g(7) - состояние упрочнения;
при у = 1/2 G + Г/2т, - имеем упругопластическую среду;
при v =
- вязкопластическую среду и т.д.
1.8. Уравнение теплопроводности.
Начальные и граничные условия
В процессах пластической де
рмации вся энергия деформации
превращается в тепло, поэтому при анализе следует учитывать изме-
нение температуры.
Дифференциальное уравнение теплопроводности для движущейся
среды с учетом тепла, выделяющегося в процессе пластической де-
формации, имеет вид [2]:
<29 _ (а2° д2& । g2eL
dt ^cbc2 ду2 dz2 J pc
(1.112)
где a = X/pc - коэффициент температуропроводности; p - плотность;
с - теплоемкость среды.
Это уравнение математически описывает перенос тепла внутри
тела, устанавливая связь между временными и пространственными
dQ д& 30
изменениями температуры 9(х, 0- Поскольку —— = —-+——о,, а на-
ar dt дх,
пряжения и скорости деформации заранее не известны, уравнение
(1.112) не может в общем случае рассматриваться без анализа движе-
ния сплошной среды и одновременного определения поля скоростей
и напряжений. Таким образом, краевая задача теории теплопровод-
ности является связанной, поскольку она объединяет теорию тепло-
проводности и механику сплошных сред. Методы решения таких за-
дач мы будем развивать дальше, а пока предположим, что из каких-
то соображений скорости и напряжения нам известны. При этом
уравнение (1.112) представляет собой линейное дифференциальное
35
уравнение в частных производных второго порядка с переменными
коэффициентами v,( х, t) при производных dQ!dxr Это уравнение от-
носится к параболическому типу и имеет в общем случае бесконечное
множество решений. Для того чтобы из этого множества выбрать
конкретное, однозначно характеризующее рассматриваемый про-
цесс, необходимо к основному уравнению присоединить краевые ус-
ловия. К ним относится начальное условие
8(х,0) = <р(х)
(1.113)
характеризующее распределение температуры внутри тела D в на-
чальный момент времени, и закон взаимодействия между окружаю-
щей средой и поверхностью тела S (граничные условия).
Рассмотрим некоторые типичные граничные условия.
Граничные условия первого рода задают распределение температу-
ры по поверхности тела в любой момент времени:
0ls=y(M,4 (1П4)
где М - точка поверхности S.
Температура у(Л/, г) может быть постоянна или может меняться
с переходом от одной точки поверхности к другой и с изменением
времени t.
Граничное условие второго рода задает плотность теплового пото-
ка в каждой точке поверхности тела как функцию времени;
<Ц$ = (1И5)
где q„ — проекция вектора теплового потока q на направление внеш-
ней нормали п к поверхности тела. Поэтому величина qn положи-
тельна, если поверхность тела охлаждается, и отрицательна - при на-
гревании поверхности.
Из закона теплопроводности Фурье q - -к grad 0 следует, что
qn =q-n- -кgrad 0 л = —к—~.
Это позволяет записать условие (1.115) в виде
(1.116)
-к— =v(A/,4
В частности, тепловой изоляции поверхности тела соответствует
условие
36
= 0.
(1.118)
on
При нагреве металла в высокотемпературных печах процесс пере-
дачи тепла в основном происходит излучением. Тепловой поток, по-
лучаемый поверхностью тела от нагретых стен и свода печи, прямо
пропорционален разности четвертых степеней абсолютных темпера-
тур поверхности излучения Те и поверхности тела Т3 (закон Стефана-
Больцмана):
-д„=рС(тс4-Г;), (1.119)
где Р - постоянная Стефана-Больцмана; с - постоянный коэффици-
ент, учитывающий условия теплообмена. Минус в левой части урав-
нения (1.119) связан с тем, что поверхность тела нагревается и qn < 0.
Если температура поверхности тела Ts значительно меньше тем-
пературы печи Т, = const, то вычитаемым можно пренебречь и запи-
сать следующее граничное условие:
?я Is =РсТс4= const. (1.120)
Граничное условие третьего рода описывает процесс конвективно-
го теплообмена между поверхностью тела и окружающей средой по
закону Ньютона:
Чп Is =а(е, “О
где в, - температура окружающей среды; а - коэффициент пропор-
циональности, называемый коэффициентом теплообмена, Вт/м • град.
Эта величина численно равна количеству тепла, отдаваемого (или
получаемого) единицей площади поверхности тела в единицу време-
ни при разности температур между поверхностью и окружающей
средой в 1°С.
Воспользуемся зависимостью (1.116), что позволяет записать это
соотношение следующим образом:
-А:^ = а(0,-еД (1.121)
on
Таким образом, граничные условия третьего рода описывают
процесс теплопередачи, при котором количество тепла, передаваемо-
го в единицу времени с е;
HI
ады поверхности тела в окружающую
среду, прямо пропорционально разности между температурой по-
верхности тела и температурой окружающей среды. Однако во мно-
37
гих случаях этот закон теплопередачи применим и для описания на-
гревания или охлаждения тел лучеиспусканием. Для того чтобы убе-
диться в этом, запишем правую часть соотношения (1.119) (закон
Стефана-Больцмана) в следующем виде
рс(гс2-г;)=ре(тс + TS )(г2-Т/)=-а(Т)(0,-ес), (1.122)
где коэффициент лучистого теплообмена а(7), имеющий ту же раз-
мерность, что и коэффициент конвективного теплообмена, зависит
от температуры, свойств поверхности тел, участвующих в лучистом
теплообмене, и равен
а(г)=₽с(т; +п)(т2-т2).
(1.123)
Во многих задачах температура Тс = const и величина Ts изменяет-
ся в сравнительно узких пределах, это позволяет приближенно при-
нять коэффициент а(7) постоянным, т.е. перейти к граничному усло-
вию третьего рода.
Граничное условие четвертого рода описывает теплообмен по-
верхности тела с омывающей его жидкостью или газом, а также теп-
лообмен соприкасающихся твердых тел в предположении, что темпе-
ратура соприкасающихся поверхностей одинакова.
При этом наряду с равенством температур
0, =
(1.124)
имеет место и равенство тепловых потоков:
(1.125)
Очевидно, граничное условие четвертого рода более точно, чем
условие третьего рода, описывает процесс контактного теплообмена,
особенно для нестационарных температурных полей. Вместе с тем,
применение этого граничного условия предусматривает решение еще
одной температурной задачи - анализ распределения температуры
в окружающей среде. Поскольку значения теплофизических констант
поверхностного теплообмена для реальных процессов, когда поверх-
ность тела покрыта слоем окислов или смазки, известны нам лишь
приближенно, во многих случаях приходится ограничиваться при-
ближенным описанием условий на поверхности, граничными уравне-
ниями первого или третьего рода. Рассмотрим в связи с этим сле-
дующие задачи.
38
1. Тепловой контакт нагретой заготовки и инструментов
Исследуя условия теплообмена вблизи контактной поверхности
между инструментом и нагретой заготовкой, можно сформулировать
следующую задачу. Два полуограниченных тела D\ и D2, равномерно
нагретых до температуры соответственно О1о и 02О» при t = 0 приве-
дены в соприкосновение. Найти распределение температуры в любой
момент времени.
Решение. Расположим систему координат, как указано на
рис. 1.8. При этом начало координат находится на поверхности кон-
такта. Индекс 1 при обозначении переменных будем относить к телу Df,
а индекс 2 - к телу Z>2.
Рис. 1.8. Идеальный тепловой контакт
2
х2, с2,
Температура в любой точке рассматриваемых тел зависит только
от двух переменных: координаты х и времени t. Запишем уравнения
теплопроводности:
00! (х, z) 020.(х, z) /
для тела = %,----(t > 0; х > 0);
ст дх
_ 002(х, z) 0202 (х, t) /
для тела D, —< - т,----------(Z > 0; х < О
dt 2 дх2
рмулируем краевые условия. К ним относятся начальные условия
0j (х, 0) = 01О, 02 (х, 0) = 02О и граничные условия четвертого рода:
01(+O,z) = 02(-O,/X
301(0,0 к2 ае2(о,г)
дх к^ дх ’
(1.126)
(1.127)
где 0j(+O,z) и 02(-O,z) означают предельные значения температуры
при стремлении х к нулю соответственно справа или слева.
39
Добавим дополнительное условие
описывающее затухание тепловых возмущений на бесконечности.
Решение поставленной краевой задачи с применением интеграль-
ного преобразования Лапласа имеет следующий вид:
где Дб = 01О ~е20.
Коэффициент К характеризует «тепловую активность» первого
тела по отношению ко второму. Он равен:
(1.130)
где к — коэффициент теплопроводности; с - теплоемкость; р - плот-
ность среды; b = -Jkcp - коэффициент тепловой активности.
Индексы 1 и 2 относятся соответственно к телу Di или D2-
Функция erf(z) представляет собой интеграл вероятности:
du\ erf(O) = O.
(1.131)
Таблицы этой функции имеются в математических справочниках
и курсах теории вероятности. Наконец,
erfc(z) = 1 -erf(z); erfc(O) = 1. (1.132)
Рассмотренная задача может быть использована для приближен-
ного анализа распределения температуры в окрестности контакта
между нагретым телом и инструментом в процессах обработки ме-
таллов давлением. Важным следствием из полученного решения яв-
ляется то обстоятельство, что на поверхности контакта температу-
ра сразу же после соприкосновения устанавливается равной
40
0j(+O,r)- 02(~O>O-^2O + , . ~
1 + A.
(1.133)
и остается постоянной на протяжении всего процесса теплообмена.
Это позволяет свести граничные условия четвертого рода к гранич-
ным условиям первого рода, существенно упрощая решение задачи
теплообмена.
Если тепловая активность первого тела, т.е. величина
значительно больше тепловой активности другого, то К » 1 и темпе-
ратура на контакте близка к начальной температуре первого тела.
В противном случае К —> 0 и по формуле (1.133) температура на по-
верхности соприкосновения близка к начальной температуре второго
тела. Если тепловые активности одинаковы, К = 1 и температура на кон-
такте равна среднему арифметическому начальных температур в10 и 62о-
Устремив в выражениях (1.128) и (1.129) время t к бесконечности,
найдем, что в стационарном состоянии температура во всем объеме
обоих тел будет одинаковой, равной температуре на контактной по-
верхности.
2. Тепловой контакт между инструментом и нагретой заготовкой
при наличии промежуточного слоя окалины или смазки
Предыдущая задача описывает тепловое взаимодействие между
нагретой заготовкой и инструментом для случая идеального тепло-
вого контакта, который имеет место лишь в сравнительно редких
случаях (например, прокатки заготовок с неокисленной поверхно-
стью в вакууме). Для того чтобы оценить влияние слоя окалины или
смазки между контактирующими телами, несколько изменим условия
этой задачи, предположив, что на плоскости х = 0 имеется сопротив-
ление потоку тепла между соприкасающимися телами.
Решение. Сохраняя принятую ранее систему обозначений (рис. 1.9),
заменим при постановке краевой задачи уравнения (1.126) и (1.127) ус-
ловием теплового сопротивления
Рис. 1.9. Тепловой контакт
при наличии промежуточного слоя
где коэффициент at2 характеризу-
ет условия теплообмена при нали-
чии промежуточного слоя вещест-
ва, затрудняющего перенос тела
от поверхности первого тела к по-
верхности второго.
41
Распределение температуры для первого тела
где
/(х,й,х,0=егГс
61 = 010 ~ кГк
А + А»
X
~ ехр(йх+й2х*) erfc
(1.135)
(1.136)
Индекс 1 относится к первому телу, а приведенный коэффициент
теплообмена hi вычисляется по формуле
ai?fe + Аа)
^1^2
(1.137)
Сопоставим это решение с решением краевой задачи для получе-
ния бесконечного тела, нагретого в начальный момент времени I = О
до температуры в0 = const и охлаждаемого на поверхности х ~ 0 по
закону Ньютона (граничное условие III рода):
к— = -a(e-eJ
дх v
(1.138)
Для этого случая имеем:
в = Оо-Д9/(х,й,х,г),
(1.139)
где Дв = 0О - вс, ah = aJk - относительный коэффициент теплообмена.
Отсюда следует вывод, что решение сравнительно сложной задачи
о тепловом контакте двух полуограниченных тел при наличии тепло-
вого сопротивления между ними можно свести к решению более про-
стой задачи охлаждения одного из этих тел (граничные условия
III рода), определив относительный коэффициент теплообмена h по
формуле (1.137)
й = (1.140)
и положив температуру окружающей среды равной
®ioA +
Й)+Й2
(1.141)
42
3. Распределение тепла трения между инструментом и заготовкой
В процессах пластической деформации на контактной поверхно-
сти выделяется тепло, пропорциональное мощности, которую разви-
вают силы трения. Распределение этого теплового потока между
деформируемым металлом и инструментом заранее неизвестно. При-
ближенно сведем эту задачу к следующей. Два полуограниченных те-
ла с температурой 90 соприкасаются по плоскости х = 0. При t = 0 на
контактной поверхности начали действовать равномерно распреде-
ленные тепловые источники. Интенсивность источников, т.е. количе-
ство тепла, выделяемое на единице площади за единицу времени,
а — const. Рассчитать температурное поле в каждом теле и распреде-
ление теплового потока источников между телами.
Решение. Распределение температуры находится решением
системы уравнений:
(t>0, х>0);
^2. = X2^.(t>0, х<0)
dt 2 8х2 V
при следующих краевых условиях:
z=o ел(х,о)=е2(х,о)=е0;
(1.142)
/ л х , ч 50. 502
0^+0,/) = 02 (-0,г), ---+ fc2---= —а.
дх дх
Решение этой задачи можно представить следующим образом:
в» =0О +^7vierfc-^=(x>0); (1.143)
02 = 00 + feL (х < о), (1.144)
27x2'
где
<30
ierfcu= ferfcXdX = ~L^e~u -uerfcu, (1.145)
43
a aj и а2 - тепловой поток из источников соответственно в первое
и второе тело, так что
04 +а2
(1.146)
Из граничных условий (1.142), полагая в решении (1.143), (1.144)
х = 0, получаем:
a1/Z>1=a2/fe2. (1.147)
Решая систему (1.146), (1.147), находим следующее распределение
тепловых потоков:
a.bt
(1.148)
Оно остается неизменным во времени и полностью определяется
соотношением коэффициентов тепловой активности bi (i = 1,2) кон-
тактирующих тел.
Итак, нами поставлена краевая задача теории теплопроводности:
найти распределение температуры в области D, удовлетворяющее
уравнению теплопроводности (1.112) и краевым (начальным и гранич-
ным) условиям. В курсах математической физики доказывается, что
при некоторых дополнительных предположениях решение постав-
ленной задачи существует, является единственным и устойчивым.
Другими словами, краевая задача поставлена корректно.
Контрольные вопросы
1. Чем характеризуется напряженное соответствие в точке?
2. Какие напряжения называются главными?
3. Что называется максимальным касательным напряжением?
4. Что такое инварианты тензора напряжений? Сколько их?
5. Что такое девиатор напряжений?
6. Что называется интенсивностью касательных напряжений?
7. Чем характеризуется деформированное состояние в точке?
8. Что называется главными удлинениями?
9. Что такое максимальный сдвиг?
10. Какие деформации называются малыми?
11. Что называется интенсивностью деформаций сдвига?
12. Что такое девиатор деформации?
13. Какие инварианты у девиатора деформаций?
14. Что такое скорость деформаций?
15. Что собой представляет тензор скоростей деформаций?
44
16. Какие инварианты у тензора скоростей деформаций?
17. Что называется интенсивностью скоростей д
рмаций?
18. Какой вид имеют дифференциальные уравнения равновесия в прямо-
угольной системе координат?
19. Какой вид имеет условие текучести Мизеса-Губера?
20. Каковы соотношения между напряжениями и деформациями?
21. Что такое граничные и начальные условия?
Глава 2
ТЕОРИЯ РАЗРУШЕНИЯ
И ПЛАСТИЧНОСТЬ МЕТАЛЛОВ
В главе, следуя [2], описана гипотетическая теория деформируемо-
сти металла при его обработке давлением. Под деформируемостью
понимается способность металла претерпевать пластическое формо-
изменение без разрушения (возникновения макротрещин). Матема-
тический аппарат этой теории позволяет определить допустимые
пластические деформации, не приводящие к растрескиванию. Он ос-
нован на экспериментальных данных о зависимости пластичности от
показателя напряженного состояния или на так называемых диа-
граммах пластичности. Описаны диаграммы пластичности некото-
рых сталей и цветных металлов.
2.1. Гипотеза о разрушении металлов
при пластической деформации
Общепризнанно, что пластическая деформация металлов сопро-
вождается непрерывным образованием и развитием субмикро- и мик-
ротрещин. Процесс образования трещин (очагов разрушения) связы-
вают с движением дислокаций вследствие пластической деформации
и взаимодействием полей напряжений, окружающих дислокации.
По мнению И.А. Одинга, зародыш трещины возникает в некото-
ром микро- или субмикроскопическом объеме скопления дислока-
ций, в котором упругая энергия деформации достигла некоторой
предельной величины, равной скрытой теплоте плавления. Он счита-
ет, что такая насыщенность энергией вызывает разрушение металла.
Сам по себе зародыш трещины устойчив. Однако на его остром кон-
це опять образуется линейная дислокация, которая взаимодействует
с проходящими около нее дислокациями. Это приводит к постепен-
ному разрастанию зародыша.
Дж.Дж. Гилман отмечает, что «одним из наиболее важных меха-
низмов, при помощи которых образуются трещины в твердых телах,
является локализованное пластическое течение». Он указывает на
три различных механизма трещинообразования. Дискообразный
46
сдвиг - тонкая область пластического течения. Концентрация каса-
тельных напряжений в этой области приводит к образованию тре-
щин на границах диска. Поверхность раздела пластической дефор-
мации - также тонкая область, где пластическая деформация в теле
внезапно меняется до некоторой меньшей величины. Следующий тре-
тий способ возникновения трещин - пересечение линий скольжения.
В литературе по физике металлов имеется много эксперименталь-
ного материала, подтверждающего положение о том, что пластиче-
ская деформация сопровождается возникновением и развитием тре-
щин. Это положение о возникновении и развитии трещин, на наш
взгляд, следует положить в основу теории разрушения металла при его
обработке давлением, т.е. в процессах значительного формоизменения.
Введем некоторую скалярную величину V - трещиноватость, ко-
торая будет характеризовать пораженность элементарного объема,
окружающего данную частицу металла микродефектами. Эта катего-
рия, следуя высказыванию М.М. Филоненко-Бородича, будет основ-
ной рабочей гипотезой расчетного аппарата, который построим ни-
же. Будем считать, что приращение V, полученное за некоторый
промежуток времени, пропорционально приращению степени пла-
стической деформации и обратно пропорционально пластичности Л.р
металла при реализуемом в данный момент напряженном состоянии,
показатель которого к =с!Т. Таким образом,
dA
9
или
(2.1)
Использование такой концепции при построении теории разру-
шения не ново. Понятие о поврежденное™ материала известно из
теории усталостного разрушения. Деформация тела знакоперемен-
ными напряжениями, не превышающими предел упругости, приводит
к накоплению и развитию повреждений - микротрещин. Каждый
цикл знакопеременного напряжения образует некоторую долю по-
врежденности Д£>,. После действия п циклов поврежденность
Л
D„ = У ДГ>.
1
Считают, что к моменту образования трещины усталости (и = N)
поврежденность равна единице (DN =1).
47
Аналогичное понятие известно и в теории ползучести. Л.М. Кача-
нов, рассматривая вопрос о времени до разрушения в условиях пол-
зучести, ввел понятие о некотором скаляре поврежденности 1|/. В на-
чальный момент, когда еще нет заметных пластических деформаций,
автор условно принял у = 1. В момент исчерпания пластичности
и хрупкого разрушения у = 0.
В теории разрушения принято, что период накопления микротре-
щин составляет основную часть времени «жизни» (деформации) об-
разца. Процесс образования макротретин усталости протекает сравни-
тельно быстро. В упомянутой выше теории разрушения при ползучести
Л.М. Качанов, ссылаясь на работу Я.Б. Фридмана, также различает две
стадии протекания разрушения: постепенное развитие трещин и значи-
тельное ускорение разрушения, причем преобладает первая стадия.
Также на две стадии можно разделить процесс разрушения при
больших пластических деформациях. По мере развития деформации
растут зародыши трещин, начинают все сильнее действовать эффек-
ты концентрации напряжений. До некоторых пор трещина остается
устойчивой, и для ее дальнейшего развития необходимы дополни-
тельные пластические деформации. Вторая стадия начинается с неко-
торого момента, когда трещина достигает критического размера
и теряет устойчивость. Вопросы самопроизвольного развития тре-
щин рассматривались Гриффитцем, Орованом, Ирвином, Баренблат-
том, Качановым, Друккером и др. После достижения критического
размера достаточно небольшой пластической деформации, чтобы
трещина резко увеличила свои размеры. Несколько трещин объеди-
няются, образуя поверхность разрушения. Степени деформации эле-
ментарного объема в начале и в конце второй стадии отличаются не-
значительно, и теория разрушения, на наш взгляд, может рассматри-
вать лишь первую стадию.
Параллельно с процессом возникновения и увеличения микроде-
фектов (трещин) в пластически деформируемом теле идут процессы
«залечивания» зачатков нарушения сплошности и торможения их
развития. Соприкосновение поверхностей трещины в условиях сжа-
тия и их относительное перемещение из-за пластической деформации
могут вызвать схватывание (сварку).
Явление залечивания хорошо известно по работам сотрудников
Института физики металлов АН СССР. Путем пластического растя-
жения образцов из чистой меди они получили металл, структура ко-
торого была поражена микротрещинами и микропорами. Пластич-
ность его была в значительной степени исчерпана. Затем эти же об-
разцы были продеформированы в условиях всестороннего сжатия.
48
При этом наблюдалось полное залечивание пор и трещин, а пластич-
ность повысилась до первоначального уровня. Авторы указали, что
всестороннее сжатие образцов без их пластической деформации не
привело к сколько-нибудь заметному уменьшению количества и раз-
меров дефектов.
Итак, уменьшение пораженности элементарного объема микро-
дефектами за промежуток времени dx можно записать:
=-C2^dt.
(2.2)
Коэффициент Q зависит от напряженного состояния. Если с/Г > 0, то
с2 = 0. Залечивание происходит лишь при о/Г < 0.
Поскольку процессы развития микродефектов и их залечивания
могут идти одновременно, то результирующее приращение
Hdx
(2.3)
Соотношение неизвестных пока величин ci и с2 определяет интен-
сивность охрупчивания (dy > 0) или улучшение пластических свойств
(<А|/ < 0) в результате малого акта пластической деформации со сте-
пенью деформации сдвига Hdx.
Моменту разрушения соответствует некоторая трещиноватость Ур.
Поделим правую и левую часть соотношения (2.3) на и обозначим
выражение в круглых скобках через В, а У/у р = V - степень исполь-
зования запаса пластичности. Получим
Hdx
(2-4)
Предположим, что правая часть не зависит от у . После интегри-
рования дифференциального уравнения (2.4) имеем
Hdx
(2.5)
Формула (2.5) показывает степень использования запаса пластич-
ности металла к моменту времени Л Если иметь в виду, что в исход-
ном недеформированном состоянии у = 0, а в момент разрушения
49
V = 1, то условие де
пТк
рмирования за период времени 0...Г элементар-
ного объема металла без разрушения имеет вид
Н(т)
(h<l
(2.6)
где В, Н, Лр, к, - известные функции, описывающие напряженно-
деформированное состояние некоторого элементарного объема ме-
талла и его пластические свойства.
Условие деформирования без разрушения (2.6) не учитывает зале-
чивание микродефектов за счет явлений рекристаллизации и за счет
диффузионных процессов, которые происходят при высоких темпе-
ратурах. Поэтому условие (2.6) далее будет применяться для оценки
предельных деформаций металлов в условиях холодной обработки
давлением.
Распространим описанную выше теорию разрушения на процессы
горячего деформирования. Как уже отмечалось, высокая температу-
ра горячей обработки давлением способствует залечиванию дефек-
тов, возникающих при пластической де
jus
рмации. Это залечивание
идет во времени. Чем больше времени прошло после деформации,
тем полнее восстанавливается пластичность. В отличие от холодной
деформации металл уже не обладает «идеальной памятью». Чтобы
учесть это обстоятельство, в условие деформирования без разруше-
ния (2.6) под интеграл следует ввести коэффициент наследственности
E(t - т), который изменяется от 0 до 1 и является монотонно убы-
вающей функцией аргумента.
Таким образом, для процессов горячей обработки металлов дав-
лением условие деформируемости можно принять
у = f Е(/ - т)в(т) -т-гт-
J ЛД£(т)]
rfxcl.
(2.7)
При холодном деформировании коэффициент наследственности
Е - 1. Проследим за деформацией некоторого элементарного объема
металла. Пусть формоизменение происходит при постоянном пока-
зателе напряженного состояния и процесс деформации близок к мо-
нотонному. Тогда, вероятно, между скоростью развития трещин
и скоростью их залечивания установится неизменное соотношение
(В = const) и условие (2.7) примет вид
50
— [ял<1 или ЯЛ <Л
Л J р
р о
Если теперь примем, что В = 1, то придем к очевидному условию
деформируемости: степень деформации в условиях формоизменения,
близких к монотонному, с постоянным о/ Т не должна превышать не-
которую предельную величину, свойственную данному напряженно-
му состоянию. Таким образом, в случае деформирования некоторого
элементарного объема в условиях формоизменения, близких к моно-
тонному процессу, с постоянным показателем напряженного состоя-
ния В - 1.
Одной из основных величин, входящих в условия деформирова-
ния без разрушения (2.6) и (2.7), является Лр, точнее, зависимость Лр
от показателя напряженного состояния <з!Тъ условиях, близких к мо-
нотонному деформированию, т.е. диаграмма пластичности, которая
должна быть определена экспериментальным путем.
Контрольные вопросы
1. Какие понятия использует теория накопления повреждений?
2. Как записывается условие деформируемости для холодной дефор-
мации?
3. Как записывается условие деформируемости для горячей деформации?
51
Глава 3
МАТЕМАТИЧЕСКОЕ
МОДЕЛИРОВАНИЕ ПРОЦЕССОВ
ОБРАБОТКИ МЕТАЛЛОВ
ДАВЛЕНИЕМ
Математическая модель процесса обработки металлов давлением
представляет собой систему уравнений, позволяющую найти в облас-
ти пластического течения распределение скоростей, деформаций,
температур, напряжений, рассчитать вероятность разрушения метал-
ла и определить оптимальные условия его деформирования.
Для получения единственного решения, соответствующего рас-
сматриваемому технологическому процессу, необходимо задать на-
чальные и граничные (краевые) условия, т.е. поставить краевую за-
дачу обработки металлов давлением.
3.1. Краевая задача обработки металлов давлением
Краевая задача обработки металлов давлением включает замкну-
тую систему уравнений, описывающую движение сплошной среды
в области D с границей 5, и соответствующие начальные и граничные
условия.
Система уравнений движения сплошной среды, приведенная ра-
нее, может быть выражена следующим образом.
Деформационная теория:
уравнения равновесия:
a.7J=0; (3.1)
уравнение связи деформаций с перемещением:
(3.2)
уравнение связи тензора напряжений с тензором деформаций:
Г = 2уТ,
где
52
Теория течения-.
уравнения равновесия ст... = 0;
уравнения связи скорости деформаций со скоростью частиц
(3-4)
уравнение связи тензора напряжений с тензором скоростей де-
формаций
Я = 2уТ,
(3.5)
Значения функции \|/ зависят от реологической модели среды. Для
ряда моделей они были определены выше.
Начальные и граничные условия. Поставить конкретную задачу
о движении сплошной среды в области D с границей S означает вы-
брать математическую модель среды, т.е. записать соответствую-
щую замкнутую систему уравнений, задать внешние силы и сформули-
ровать начальные и граничные условия.
К начальным условиям относятся уравнения, описывающие рас-
пределение искомых давлений, температур, скоростей в начальный
момент времени. В некоторых случаях одних только начальных ус-
ловий вполне достаточно для выделения определенного решения
(например течение в неограниченной области).
Примером задания начальных условий является широко исполь-
зуемая гипотеза о первоначальном напряженном состоянии тела.
Условия на границе S (известной или неизвестной при любом t)
можно разделить на механические и температурные. К последним
относятся уравнения, описывающие распределение температуры или
условия теплоотдачи на границе тела.
Температурные граничные условия. При задании температурных
граничных условий наибольшие трудности связаны с описанием
процесса теплопередачи между пластически деформируемым телом
и инструментом. Строгий подход к этой проблеме предусматривает
совместный расчет температурных полей для этих тел, т.е. решение
двух достаточно сложных задач с согласованием решений на границе
раздела.
Вместе с тем, во многих случаях, учитывая отсутствие полной ин-
формации о значениях теплофизических констант на контактных по-
53
верхностях, а также сложное влияние на процессы теплопередачи
промежуточных пленок окислов или смазок, можно ограничиться
приближенным описанием граничных условий, воспользовавшись
результатами решения сравнительно простых задач о контакте полу-
ограниченных тел. При этом методика задания граничных условий
сводится к следующему. Вводится усредненная по объему тел на-
чальная температура обрабатываемого металла 0ю и инструмента 0го
(здесь и в дальнейшем индексы «1» и «2» относятся соответственно
к обрабатываемому металлу и инструменту). Теплообмен на кон-
тактной поверхности моделируется теплообменом двух полуограни-
ченных тел.
Наиболее важными теплофизическими константами, определяю-
щими условия теплообмена, являются коэффициенты тепловой ак-
тивности металла и инструмента 62 = д/^г^гРг > где _
коэффициент теплопроводности; с - теплоемкость; р - плотность
среды.
Рассмотрим три наиболее важных варианта граничных условий.
1. Идеальный тепловой контакт между металлом и инструмен-
том. Граничные условия (ГУ) сводятся к ГУ I рода:
e|s=ea>+#F- <3.6)
1 4- А
где де = Ою...020, К = 6i/62 •
2. Теплообмен при наличии теплового сопротивления между кон-
тактными поверхностями, связанного с наличием промежуточного
слоя окалины или смазки. Граничные условия записываются как ГУ
II рода:
где относительный коэффициент теплообмена равен
а(^+г>2)
(3-8)
а температуру окружающей среды определяет формула
6»Д +
(3-9)
54
3. Теплообмен при наличии теплового сопротивления и выделения
тепла трения между контактными поверхностями Граничные усло-
вия сводятся к ГУ III рода:
= й(0-0'),
(3.10)
где h вычисляется по формуле (3.8), а
о'=ес
оД
+ й(^+Л2)’
(3.11)
где 0С определено зависимостью (3.9); а- интенсивность тепловых
источников на границе тела.
Механические граничные условия. Механические граничные усло-
вия могут быть динамическими (статическими), кинематическими
и смешанными.
Динамические граничные условия задают на границе тела векто-
ром напряжений б" как известную функцию точки границы М
и времени V.
ал = (3.12)
Так, при движении сплошной среды можно рассматривать по-
верхности, называемые свободными, на которых поверхностные на-
пряжения сводятся к атмосферному давлению.
Если среда находится в равновесии, то условия (3.12) называются
статическими.
Если среда является жесткопластической, то динамические
(статические) условия на границе между жесткой областью и обла-
стью пластического течения могут задавать выраженные через по-
верхностные напряжения главный вектор Q и главный момент М
относительно некоторой точки О известных согласно условиям зада-
чи внешних сил, приложенных к жесткой области (рис. 3.1, а).
Кинематические граничные условия полностью определяют на по-
верхности тела вектор перемещения и или скорости и.
Пусть нам задано положение и движение некоторого участка гра-
ницы. При отсутствии проскальзывания материальных частиц по ка-
сательной к границе имеет место условие «.прилипания» частиц к по-
верхности. Это условие запишется в виде:
(313)
^ср ^пов'
55
Рис. 3.1. Условия на границе между жесткой областью
и областью пластического течения
Граничное условие прилипания к обтекаемым телам характерно
для вязкой жидкости. Эксперименты показывают, что в определен-
ных условиях (например при прокатке высоких полос) прилипание
имеет место на контакте между инструментом и горячим металлом.
К этому же условию сводится условие непрерывной стыковки дви-
жения «жесткой» области с полем скоростей в области пластического
течения для жесткопластических сред (рис. 3.1, б). В этом случае мо-
гут быть заданы скорости поступательного движения и и углового
вращения «> жесткой зоны.
Если частицы среды могут скользить вдоль границы, то на по-
верхности действует более слабое кинематическое условие обтекания
(непроницаемости):
^ЛСр ^ЛПОВ’ ^лср ^ЯПОВ* О-14)
При этом на границе совпадают только нормальные составляю-
щие скорости (перемещения) частиц среды и граничной поверхности.
Это условие дополняется некоторым физическим законом внешнего
трения на контактной поверхности. Так, для идеальной жидкости ка-
сательные напряжения трения равны нулю; для вязкой жидкости они на-
столько велики, что приводят к прилипанию частиц к поверхности.
В общем случае закономерности внешнего трения накладывают
некоторые условия на соотношения между нормальной и касатель-
ной составляющими вектора поверхностных напряжений.
Таким образом, мы пришли к смешанным граничным условиям,
в которых на поверхности тела частично накладываются ограничения
на кинематику и частично на поверхностные напряжения. Остановим-
ся более подробно на очень важном вопросе трения на контактной по-
верхности.
56
Внешнее трение при обработке металлов давлением. Как известно,
внешнее трение - это механическое сопротивление, возникающее на
поверхности касания двух соприкасающихся тел при их относитель-
ном перемещении. Различают два вида трения: качения и скольжения.
Для обработки металлов давлением характерно трение скольжения,
которое по своей природе существенно отличается от трения сколь-
жения в узлах машин. Отсылая читателя, желающего ознакомиться
более детально с этим вопросом, к рекомендуемой литературе, ко-
ротко сформулируем основные результаты многочисленных экспе-
риментальных исследований.
Рассмотрим элементарную площадку, расположенную на кон-
тактной поверхности, с нормалью Я. Вектор поверхностных напря-
жений а", действующий на площадке, представим в виде суммы двух
векторов: вектора нормального давления р и вектора напряжения
трения т (рис. 3.2):
ол = Я + т. (3.15)
Вектор р направлен по внут-
ренней нормали к поверхности,
его модуль р равен нормальной
силе, действующей на единицу
площади.
Вектор т лежит в плоскости, ка-
сательной к поверхности, и направ-
лен в сторону, противоположную
вектору скольжения частиц металла
относительно инструмента. Пусть
и - скорость частицы, б] - ско-
Рис. 3.2. Нормальные
и касательные напряжения
на контактной поверхности
рость инструмента на рассматри-
ваемой площадке, Д и = и - ut - скорость скольжения. Тогда
Ди
]дб|’
(ЗЛ6)
где т - модуль вектора напряжения трения; обычно называется про-
сто напряжением трения.
Эксперименты свидетельствуют о сложном характере зависимости
напряжения трения от основных факторов: давления р, скорости
скольжения До, суммарного относительного перемещения Дй, со-
стояния поверхности инструмента, химического состава и темпера-
57
туры инструмента и деформируемого тела, наличия смазки или ока-
лины на контактной поверхности. Практически используются два
упрощенных закона трения:
1) закон Амонтона-Кулона
(3.17)
где/- коэффициент трения;
2) закон Прандтля (Зибеля)
т =А,
(3.18)
где а, - предел текучести обрабатываемого металла, а коэффициент f
(его тоже часто условно называют коэффициентом трения), как
и в формуле (3.17), должен отражать влияние перечисленных выше
факторов.
Уравнение (3.17) обычно используется для описания внешнего
трения при холодной обработке металлов давлением, а уравнение
(3.18) - при горячей обработке металлов.
Необходимо отметить, что в соответствии с условием текучести
напряжение трения не может превышать предела текучести сдвига т,
деформируемого металла в приконтактном слое. Поэтому уравнение
(3.17) более правильно записать следующим образом:
fp при fp<y,
Ь при fP^ts
(3.19)
Вследствие интенсивной пластической деформации приконтакт-
ных слоев предел текучести в этой области может существенно пре-
вышать усредненное по всей области пластического течения значение
этой величины.
Краевые задачи механики сплошных сред обладают специфиче-
скими особенностями. К ним относится нелинейность основных
уравнений, сложная геометрия области течения, наличие заранее не-
известных участков границы, определенные трудности при описании
граничных условий. Это приводит к стремлению упростить поста-
новку задачи и без внесения существенных погрешностей описать
лишь наиболее важные стороны исследуемых явлений.
Типичные упрощения. К упрощениям, используемым в постановке
краевых задач механики сплошных сред, относятся уменьшение чис-
ла независимых переменных и методы линеаризации.
В некоторых практически приемлемых случаях в соответствую-
щей системе координат движение сплошной среды и протекающие
58
в ней физические процессы можно считать установившимися (ста-
ционарными). Используя переменные Эйлера, можно исключить вре-
мя г, т.е. сократить число независимых переменных на единицу. Су-
щественно упрощает решение отсутствие в таких задачах начальных
условий, поскольку во всех уравнениях для стационарных процессов
выпадают частные производные по времени.
Во многих задачах при правильном выборе системы декартовых
координат движение сплошной среды можно считать плоским. Ме-
тоды решения соответствующих «плоских» задач гидромеханики,
теории упругости и пластичности хорошо разработаны. В этом слу-
чае удается привлечь мощный аппарат и методы теории функций
комплексного переменного, метод характеристик.
В задачах с наличием осевой симметрии использование цилиндри-
ческой системы координат позволяет исключить из рассмотрения уг-
ловую координату о, оставив вещественными аргументами только
координаты г, z и t. При этом все уравнения и формулы, дающие ре-
шение, будут инвариантны относительно поворотов на любой угол
вокруг оси z.
Методы линеаризации позволяют свести решение нелинейных за-
дач к итерационному процессу, в котором на каждом этапе решается
некоторая линейная задача.
3.2. Методы решения краевых задач
Разрабатывая математическую модель процесса, происходящего
в сплошной среде, например пластического течения металла при об-
работке давлением, после постановки краевой задачи необходимо
перейти к математическому исследованию процесса. Для этого нужно
получить решение краевой задачи — точное или приближенное.
Как правило, точное решение удается получить при самой упро-
щенной постановке задачи для наиболее грубых и несложных моде-
лей. На практике приходится иметь дело с течением среды в области
сложной геометрии с неравномерным распределением температуры.
Большие трудности связаны также с описанием реологических
свойств среды, моделирующей реальные свойства металлов и спла-
вов, условий трения и теплопередачи на контактной поверхности те-
ла. Все это приводит к необходимости применения приближенных
методов и решения задач на современных электронных вычисли-
тельных машинах (ЭВМ), выполняющих сотни тысяч и даже миллио-
ны операций в секунду.
59
Наиболее перспективными методами, решающими указанную
проблему, являются методы конечных и граничных элементов. Они
и будут рассмотрены ниже.
Контрольные вопросы
1. Что включает в себя математическая модель процессов ОМД?
2. Как формулируется краевая задача ОМД по деформационной теории?
3. Как формулируется краевая задача ОМД по теории течения?
4. Как формулируются начальные и граничные условия?
5. Какой вид имеют температурные граничные условия?
6. Что включают в себя механические граничные условия?
Глава 4
МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ
ДЛЯ РЕШЕНИЯ ЗАДАЧ
ТЕПЛОПРОВОДНОСТИ
4.1. Основные концепции метода конечных
элементов (МКЭ)
В этой главе, следуя [3], рассмотрим сущность МКЭ на примере
решения двумерной стационарной задачи теплопроводности в облас-
ти D произвольной формы (рис. 4.1):
д
дх
дТ'
дх >
3x1 ду,
+?и =0;
(4.1)
при граничных условиях третьего рода:
(4.2)
где Л - теплопроводность; а - коэффициент теплоотдачи на границе;
qv, qs - объемная и поверхностная плотности мощности источников
теплоты.
Величины X, a, q^ qs могут быть заданы в виде произвольных
функций коор
IIU
чат х, у, в том числе и кусочно-непрерывных функций.
Метод конечных элементов основан на определении температур-
ного поля путем приближенного решения соответствующей вариаци-
онной задачи. Для формулировки этой задачи напомним понятие
функционала. Оператор 7[Дх)] называется функционалом, заданным
на некотором множестве функций, если каждой функции/(х) из этого
множества по некоторому правилу ставится в соответствие числовое
значение 7[/(х)]. Иными словами, функционал является как бы
«функцией от функции». В практических приложениях обычно встре-
чаются функционалы, заданные в виде некоторых интегралов, в по-
дынтегральные выражения которых входят функции Дх) .
61
Рис. 4.2. Разбиение области D на
конечные элементы
Рис. 4.1. Область для решения
двумерной задачи теплопроводности
Многие краевые дифференциальные задачи теплопроводности
и конвективного теплообмена эквивалентны задачам отыскания
функций, доставляющих минимум некоторым специально сконст-
руированным функционалам. Задача на отыскание функции, мини-
мизирующей функционал, называется вариационной. На основе
перехода от краевых дифференциальных задач к вариационным
развиты многие приближенные аналитические методы решения за-
дач теплопроводности. Отметим, что возможность вариационной
формулировки задачи определения температурного поля (4.1), (4.2)
обусловлена свойствами дифференциального оператора уравнения
теплопроводности. Мы приведем вариационную формулировку
рассматриваемой краевой дифференциальной задачи (4.1), (4.2) без
доказательства. Задача решения уравнения (4.1) с граничными ус-
ловиями (4.2) эквивалентна задаче определения функции 7\х,у),
минимизирующей функционал I [ Т(х, у)] вида
L
(4.3)
Наиболее широко распространенный прием получения прибли-
женного решения вариационных задач состоит в следующем. При-
ближение для искомой функции 1\х, у) разыскивается в виде
м
Т(х,у) №^o-mfm{x,y\
m=l
(4-4)
где ат - неизвестные постоянные коэффициенты, а Д(х, у) - извест
ные функции пространственных координат. Если подставить (4.4)
62
в функционал (4.3), то можно провести интегрирование по простран-
ственным переменным и получить величину I, зависящую уже не от
неизвестной функции, а от неизвестных коэффициентов разложения
(4.4):
/ = /(«!, (4.5)
Очевидно, что для определения приближенного решения вариа-
ционной задачи в виде (4.4) следует найти значения oq, обес-
печивающие минимум функции (4.5), т.е. задача сводится к отыска-
нию точки минимума «обычной» функции нескольких переменных.
Как известно, значения eq,а м, обеспечивающие минимум функ-
ции /, находятся из решения системы уравнений
д! д!
= 0, ...,-= 0. (4.6)
ба,---------------------------------дам
Решив систему (4.6), можно найти значения а,, ...,ам, и, подста-
вив их в (4.4), определить приближенное решение вариационной за-
дачи.
Центральным местом в изложенном методе является назначение
координатных функций разложения (4.4) fi, Метод конечных
элементов основан на использовании описанной схемы приближен-
ного решения при специфическом выборе вида координатных функ-
ций /ь Благодаря этому выбору неизвестные коэффициенты
в разложении (4.4) приобретают ясный физический смысл.
Построение координатных функций проводится в МКЭ после
разбиения области определения искомой непрерывной величины на
N подобластей, называемых элементами, и фиксации в них М узло-
вых точек, выбираемых на границах элементов (например рис. 4.2).
Отметим, что число членов в разложении (4.4) равно числу узловых
точек. Каждая из функций f„(x, у) обладает следующими свойства-
ми. Значение функции /от(х, у) в т-й узловой точке с координатами
х = У ~ Ут равно единице, а в остальных узловых точках - нулю.
Кроме того, функция fm(x, у) может быть отлична от нуля только
в элементах, содержащих т~й узел. В остальной части области D она
считается равной нулю.
При таком выборе координатных функций fm(x, у) любой неиз-
вестный коэффициент ат в разложении (4.4) равен приближен-
ному значению температуры ат в т-й узловой точке. Действи-
тельно, при подстановке в аппроксимацию (4.4) координат т-го
63
узла (х = хт,у=ут) значения всех координатных функций, кроме /и-й
функции, будут равны нулю, а значение т-й функции - единице
и, следовательно:
м
rn}a Um ~ ~ (4.7)
v»l
При использовании разложения (4.4) в каждой точке области D
«работают» только те координатные функции, у которых коэффици-
енты равны приближенным значениям температур узловых точек ко-
нечного элемента, содержащего данную точку.
Отметим два важных факта, вытекающих из рассмотренных
свойств координатных функций.
Во-первых, функция l(alt получившаяся при подстановке
разложения (4.4) в функционал (4.3), будет функцией от неизвестных
приближенных значений температуры в узловых точках щ,..., им:
7 = 7(u1,...,ma/),
(4.8)
а уравнения, вытекающие из условия минимума (4.6)
д! д/
----= 0,...,---= 0, (4.9)
дим
будут алгебраическими уравнениями разностной схемы метода ко-
нечных элементов относительно искомых температур в узлах.
Во-вторых, пространственное распределение температуры внутри
любого элемента аппроксимируется суммой произведений координат-
ных функций на коэффициенты, равные приближенным значениям
температуры в узловых точках, принадлежащих данному элементу.
Координатные функции fm(x, у), т = I, М, строятся на основе так
называемых функций формы элементов. Каждая из функций формы
конкретного элемента равна единице в одной «своей» узловой точке,
принадлежащей данному элементу, и нулю в остальных узлах этого
элемента, т.е. для элемента вводится столько функций формы, сколь-
ко в нем содержится узлов. Вне элемента все его функции формы счи-
таются равными нулю. Таким образом, функция формы л-го элемен-
та, равная единице в принадлежащей ему т-й точке, является
«представителем» координатой функции fm(x, у) в этом л-м элементе.
Поэтому температурное поле в n-м элементе аппроксимируется сум-
мой произведений его функций формы на приближенные значения
температур в его узловых точках. Очевидно, что для каждого элемен-
64
та получается своя аппроксимация, но на границах элементов долж-
на сохраняться непрерывность температурного поля.
Перейдем к реализации изложенной методики для задач (4.1), (4.2).
4.2. Построение дискретной модели
и функций формы элементов
Первый этап численного решения задачи методом конечных эле-
ментов включает выбор вида элементов и способа расположения
в них узловых точек, разбиение области на элементы и размещение
узлов, а также определение функций формы. Отметим, что эти функ-
ции существенным образом зависят от вида используемых элементов
и способа расположения узлов. При решении двумерных задач при-
меняют элементы различной формы (рис. 4.3, а-г): треугольные, че-
тырехугольные, элементы с криволинейными границами, с узловыми
точками, расположенными по границам элементов. При решении
трехмерных задач элементы могут быть выбраны в виде тетраэдров,
косоугольных параллелепипедов, призм с треугольным основанием
(рис. 4.4, a-в). Ниже ограничимся рассмотрением наиболее простых
и часто применяемых линейных треугольных элементов с узлами,
расположенными в их вершинах (рис. 4.3, а). Название «линейные»
у этих элементов обусловлено тем, что при наличии на границе трех
узлов допустимо использовать только функции формы F (х, у), ли-
нейные по координатам х и у:
F (х, у) = а + Ьх + су. (4.10)
Рис. 4.3. Элементы, применяемые для решения двумерных задач
Три коэффициента а, Ь, с определяются из трех условий: в двух
вершинах значения F (х, у) равны нулю, а в одной - единице. Приме-
нение более сложных элементов позволяет использовать для аппрок-
симации температурного поля в элементе полиномы более высоких
65
степеней. Например, для треугольного элемента с шестью узлами
(рис. 4.3, б) нужно применять более сложные по сравнению с (4.10)
функции вида
F(x,y) = а + bx + су + dxy + ex2 + fy2. (4.11)
Рис. 4.4. Элементы, применяемые для решения трехмерных задач
Отметим, что при построении функций формы, кроме сформули-
рованных выше условий на их значения в узловых точках (единица
в одном узле и ноль в остальных), необходимо обеспечивать непре-
рывность аппроксимации на границах смежных элементов.
Итак, разобьем двумерную область D на N треугольных элемен-
тов и введем М узлов во всех вершинах треугольников (рис. 4.2).
Причем вершины одних треугольников не должны находиться на
сторонах других треугольников, т.е. каждый элемент должен со-
держать три узла. Присвоим сквозную нумерацию всем элементам
(л = 1, ..., N) и всем узлам (т = 1, ..., М). Номер узла в общем случае
не связан с номером элемента. Вопрос об оптимальном способе ну-
мерации узлов обсудим далее.
На рис. 4.5 отмечен и-й элемент, содержащий вершины с номера-
ми i, j, к, которые могут принимать значения от 1 до М. Эти узлы
имеют координаты хр xjt хк по оси х и ур у„ ук по оси у.
Выполненное разбиение удобно описать для дальнейшего исполь-
зования в программах следующим образом:
• задать массивы координат всех узлов {хот}^=1 и{ут}^,в по-
рядке их нумерации;
• указать для каждого л-го элемента по три номера узлов, лежа-
щих в его вершинах;
• указать для каждого л-го элемента признак, определяющий при-
надлежность его сторон внешней границе области D.
Номера узлов в вершинах и признаки принадлежности сторон
границе будем задавать в виде следующей таблицы:
66
Номер элемента Номера узлов Признак принадлежности сторон границе
~ i « А 4 п 0, или 1, или 2
Такую таблицу можно задать в виде матрицы размером Wx4, ко-
торую будем называть индексной матрицей. Номера узлов i, j, к в ин-
дексной матрице будем записывать в порядке, соответствующем об-
ходу элемента против часовой стрелки. Признак принадлежности
сторон границе принимает значение 0, если стороны элемента не
прилегают к внешней границе; значение 1, если граничной является
Сторона между узлами i и /; значение 2, если граничными являются
две стороны - между узлами i nj и между узлами i и к. Ясно, что дру-
гих ситуаций можно избежать путем соответствующего выбора узла
i, с которого начинается запись номеров в строке индексной мат-
рицы.
Предложенный способ полностью описывает выполненное раз-
биение на элементы и нумерацию узлов и не содержит избыточной
информации.
Искомыми величинами в МКЭ являются приближенные значения
температуры ит в узлах т — 1, ..., Л/. Согласно основной идее метода,
распределение температуры в каждом элементе представляется в виде
Суммы, в которую входят три функции формы элемента, умноженные
на приближенные значения температуры в его трех узловых точках.
Поэтому распределение температуры в л-м треугольном элементе
u^"\x,y) имеет вид (рис. 4.6):
и<я)(х,у) = «,7Л<п)(х,у) + М^п)(х,у)+н^п)(х,у), (4.12)
где F$n\ Ff\ F™ ~ линейные функции координат х, у вида (4.9),
равные единице в узлах i, j или к соответственно и равные нулю
в двух других узлах (рис. 4.7). Индекс л указывает, что функции фор-
мы относятся к л-му элементу.
Таким образом, для функции формы Ff>n\x,у) должны выпол-
няться равенства
F^n4xi>yi) = 1; F^Xxj, уj ) = F^(xk,yk ) = 0.
(4.13)
Аналогичные соотношения имеют место и для двух других функ-
чий F)n)(x,у), Fkn\x,y).
67
Рис. 4.5. Элемент, содержащий
вершины jt к
Рис. 4.6. Распределение температуры
в треугольном элементе
Используя условия (4.13), можно выразить коэффициенты
(т - i.j.k) через координаты узлов ij, к и определить функции формы.
Рис. 4.7. Значение функций в узлах элементов
Приведем выражения для этих коэффициентов. Здесь и далее ин-
декс элемента (п) будем опускать, помня, однако, что функции фор-
мы различны для каждого элемента. Из условий (4.13) следует:
aj ={хкУ,-Укх,У^
* b; =(Ук -л)/25>
Су — (х( — Хд. ^)2S,
(4.14)
где S - площадь треугольника, значение которой может быть вычис-
лено через координаты узлов:
5 = (х}ук -xkyj +yjXi-yt;Xi +хкУ1 -Xjyi)/2. (4.15)
68
Градиент температуры в каждом элементе имеет постоянное
значение, и производные по координатам определяются соотноше-
ниями:
бу
(4.16)
При выводе разностных уравнений МКЭ нам понадобится вычис-
лять интегралы по площади треугольника и по отдельным сторонам
от выражений, содержащих функции формы. Приведем ряд исполь-
зуемых в дальнейшем формул. Интеграл по площади элемента от
любой его функции формы равен
||4и)(х,уХх^ = 5(п)/3,
(4-17)
т =
Как видно из рис. 4.7, этот интеграл равен объему пирамиды
с треугольным основанием и высотой, равной 1.
Интегралы по стороне L„ вычисляются по формулам:
р<(х,у)<// = ^/2;
ч
|гДх,у)Л = £,у/3; (4.18)
jFi(x,y)Fj(x,y)dl = Lgj/6,
где£,у=[(х(. -xJ+U-J/)2]
1/2
- длина стороны между узлами i и J.
Данные выражения справедливы для всех функций формы Fit F„ Fk.
4.3. Система уравнений метода конечных элементов.
Локальная и глобальная матрицы
Перейдем к следующему этапу - составлению системы разностных
Уравнений МКЭ. Система уравнений для определения температур
в узлах составляется на основании условий минимума функ-
ционала (4.8). Этот функционал можно представить в виде суммы
интегралов по всем элементам, на которые разбита область D:
69
(4.19)
(4.20)
Второй интеграл в (4.20) вычисляется лишь для граничных эле-
ментов вдоль тех их сторон, которые прилегают к внешней границе L
области D.
Условия минимума функционала можно с учетом (4.19) перепи-
сать в виде
Структура системы разностных уравнении. Рассмотрим подробнее
структуру системы разностных уравнений (4.21). Возьмем т-е урав-
нение, получающееся при дифференцировании функционала (4.19) по
значению температуры в ти-м узле. Напомним, что распределение
температуры и(И) (х, у) в любом элементе зависит только от темпера-
тур ut, и., uk в узлах этого элемента. Соответственно, и значение
функционала /л) зависит только от этих температур. Поэтому
в сумме (4.19) от ит будут зависеть только /я) тех элементов, которые
включают т-й узел. Это обстоятельство позволяет формировать сис-
тему (4.21) одним из двух способов.
Первый способ основан на переборе узлов. Он заключается в вы-
делении какого-то т-го узла, в определении элементов, содержащих
этот узел, и суммировании частных производных от функционалов
этих элементов по температуре выделенного узла. Таким образом,
при получении системы (4.21) последовательно формируется первое,
второе,..., Л/-е уравнения.
Второй способ основан на переборе элементов. В этом случае вы-
деляется какой-то л-й элемент и три производные от функционала
этого элемента по его узлам заносятся в левые части соответствую-
щих трех уравнений системы (4.21), т.е. производится формирование
сразу всей системы, которое реализуется путем постепенного
«наращивания» левых частей уравнений. Этот способ получил более
широкое распространение и рассматривается ниже.
70
Очевидно, что для его реализации понадобятся выражения для
мцизводных от функционала л-го элемента /п) по температурам его
урпли ц, ик. Для их получения поменяем порядок дифференци-
крвания и интегрирования, т.е. сначала проведем дифференцирова-
подынтегрального выражения в (4.20), а затем вычислим инте-
грал по элементу от получившейся функции. При реализации этой
Процедуры будем использовать выражение (4.12) для распределения
w (х, /)» причем индекс (и) для сокращения записей опустим. Учиты-
уяя выражения (4.16) для производных, имеем
д ( ди \2
dUj \дх J
= 2{btUi +bjUj +bkuk)bi,
+ Ср*j + ^к^к ) >
---(2quu) = 2qaFi(x,y),
ди(
(4.22)
----(2qsu) = 2q3Fi(x,y),
ди,
---(aw2 )= 2a(f,w, + FjUj + Fkuk )Ft.
dui
Будем полагать, что разбиение выполнено таким образом, чтобы ве-
личины X, qv, qs, а можно было бы считать постоянными в пределах каж-
дого элемента. Если эти величины являются кусочно-непрерывными
функциями координат, то разбиение проводят так, чтобы границы эле-
ментов совпадали с линиями разрыва. Если величины являются непре-
рывными, но резко изменяющимися функциями коор,
НЯ
яат, то нужно
построить в области их сильного изменения более густую сетку.
Вычислим сначала интеграл по площади элемента, используя вы-
ражение (4.17) для интеграла от функции формы:
ди
ди
ди.
а ГэпУ
5м, v ду J
71
Если стороны элемента принадлежат границе области, то в функ-
ционале следует учесть интеграл по этим сторонам. Пусть граничной
является сторона £.. между узлами i и j . Тогда с учетом выражений
(4.18) для интегралов по сторонам от функций формы получаем:
= 2Д-(ам, /3 + ам,- /6-9j /2). (4.24)
Здесь учтено, что Fk (х, у) = 0 на стороне Li} (рис. 4.7). Если гра-
ничной является сторона Ljk, то выражение для интеграла запишется
аналогично, но вместо следует посгавить Lik, а вместо и. - ик. Если
же к границе прилегает сторона £Л, то рассматриваемый интеграл
равен нулю, так как функция формы Ft (х, у) для узла i равна нулю на
стороне £й, и, следовательно, распределение температуры на этой сто-
роне не зависит от температуры и..
Таким образом, окончательно получаем для производной от
функционала по температуре Uj следующее выражение:
а/(п)
ди.
~'KSn\bjut -ь-Ь^и: + bibkuk +c;Uj +cicjuj +cfckuk )+
* шГ J J * *
+ аД (u; I3 + Uj /б)+ &Lik («j / 3 +
-qvSn /2-qsaLik /2, (4.25)
причем слагаемые c aL,. и aLik записывают лишь в том случае, когда
соответствующие стороны принадлежат границе. Аналогичным об-
разом можно получить и выражения для производных от /И) по тем-
пературам иик.
Проанализируем теперь с учетом (4.25) структуру системы (4.21).
Видно, что производная д№1ди, представляет собой сумму произве-
дений неизвестных температур и, ик на постоянные известные ко-
эффициенты, зависящие от координат узлов и параметров задачи,
а также постоянных известных членов, не зависящих от искомых
температур. Левые части уравнений (4.21), получаемые путем сумми-
рования частных производных, имеют такую же структуру, и, следо-
вательно, приравнивая их к нулю, мы получаем линейную систему
разностных уравнений относительно неизвестных температур узло-
вых точек.
72
Локальные матрица и вектор-столбец. Для
рмирования матри-
цы линейной системы разностных уравнении удобно записать полу-
ченные выше соотношения для частных производных функционала
и-го элемента в матричном виде. Для получения матричной записи
принято использовать так называемую локальную нумерацию узлов
и соответствующих им неизвестных температур, действующую толь-
ко в рамках каждого конкретного элемента разбиения.
Остановимся на этом подробнее. Возьмем л-й треугольный эле-
мент разбиения, имеющий три узловые точки с «глобальными» но-
мерами i, j, к, и будем условно считать в рамках м-го элемента i-й
узел - первым, у-й узел - вторым, а k-й узел - третьим. Соответственно
введем локальные номера 1, 2, 3 для неизвестных температур ир ир ик
в узлах этого элемента и будем использовать следующие обозначения:
— Uj , Uj — U2 » Wjjr
(4-26)
Отметим, что соответствие между глобальными и локальными
номерами для каждого элемента разбиения задается с помощью ин-
дексной матрицы, о которой шла речь в главе 4.2.
При использовании локальной нумерации выражения для част-
ных производных от функционала д!(п' в п-м элементе можно запи-
сать следующим образом:
= gW. U(n) - <р(л), (4.27)
дй
где - матрица размером 3 х 3; , _ столбцы из трех
ди
элементов:
(4.28)
73
Из выражений вида (4.25) вытекают следующие формулы для эле-
ментов матрицы g("' и вектора ф{”>:
giin) =(tf + <f +(aLiy +aLik )/з,
= (btbj + c>c j )^(n) + a^v /6,
g^ = (b.bk +4ck frS^+aI^/6,
F(«) e»(«)
4521 ol2 >
g<2n) =(*? +4+ (aLy + a.Lik )/з,
Й? =(*A +c..cJ1S<”> +aLik /6,
(4.29)
+ а£Л)/3;
(430)
В выражениях (4.29), (4.30) для элементов матриц g(n), q/n) сла-
гаемые, содержащие множители Ltj, L,k, Lik, следует учитывать лишь
в том случае, если соответствующая сторона элемента п принадлежит
внешней границе L.
Ясно, что матрица системы линейных уравнений относительно
неизвестных температур \ит }^=1 будет формироваться на основе мат-
риц g^, а вектор-столбец свободных членов - на основе вектор-
столбцов ф'п).
Матрицу g(n’ часто называют локальной матрицей жесткости или
локальной матрицей теплопроводности, а вектор ср - локальным
вектором нагрузок или локальным вектором тепловых потоков.
Термины «жесткость» и «нагрузка» используются исторически пото-
му, что сначала МКЭ развивался применительно к задачам прочно-
стного расчета. В задачах теплопроводности в матрицы g<") входят
теплопроводности А. и коэффициенты теплоотдачи а, а в векторы
Ф(л' - свободные члены неоднородного уравнения теплопроводности
74
д граничных условий, т.е. объемные и поверхностные плотности теп-
аового потока источников теплоты. Геометрические параметры рас-
четной области учитываются коэффициентами функций
формы элемента, а также значениями Ly, L&, LjkS'-n’.
Если рассмотреть предельно крупное разбиение, при котором вся
область состоит лишь из одного элемента (N = 1), то система уравне-
ний для определения трех неизвестных температур его узлов будет
иметь вид
g(i)u(i)={p(i)>
где g(I\ локальные матрица и вектор-столбец первого и един-
ственного элемента, т.е. в случае разбиения, состоящего из одного
элемента, его локальные матрица и вектор-столбец совпадают с мат-
•Vi
ней и вектором правых частей линейной системы уравнений МКЭ.
Глобальные матрица и вектор-столбец. В реальном случае, когда
в разбиение входят N элементов, эти матрицы и вектор-столбцы, ес-
тественно, не совпадают. И в связи с использованием термина «ло-
кальный» для матрицы и вектор-столбца элемента матрица и вектор-
столбец свободных членов системы уравнений для всей области на-
зываются глобальными. Эту систему будем записывать так:
(4-31)
где G - глобальная матрица жесткости (теплопроводности) размером
М х М\ U - вектор-столбец искомых значений температур в М узлах;
Ф - глобальный вектор-столбец нагрузок (тепловых потоков).
На первый взгляд, введение дополнительной локальной нумера-
ции неизвестных в элементах разбиения и использование матричной
формы записи (4.27) представляется излишней процедурой. Однако,
как показала практика, на самом деле это позволяет сделать более
удобной процедуру формирования глобальной матрицы G и вектор-
столбца Ф при составлении программ расчета по методу конечных
элементов.
75
4.4. Формирование глобальных
матрицы и вектор-столбца.
Решение системы уравнений МКЭ
Прежде всего отметим, что процедура построения уравнений
в МКЭ имеет важную особенность по сравнению с методом конеч-
ных разностей. При построении конечно-разностной схемы мы рас-
сматривали уравнение теплового баланса для элементарного объема,
построенного около узла сетки с номером т, и сразу получали ш-е
уравнение общей системы. В случае МКЭ в т-с уравнение системы
(4.21) входит сумма производных от функционалов вычислен-
ных для различных элементов, которые содержат узел с номером т.
Поэтому при составлении каждого уравнения надо производить
суммирование «вкладов» от разных элементов. Из-за этой особенно-
сти процедура построения системы уравнений МКЭ несколько менее на-
глядна, чем в случае конечных разностей, и при ее первоначальном изу-
чении возникают некоторые трудности. Для простоты изложение начнем
с разбора конкретного примера для области, изображенной на рис. 4.8
и состоящей всего из трех элементов, которые содержат пять узлов.
Пример построения системы разностных уравнений. Нумерация
элементов, глобальная и локальная нумерация узлов приведены на
рис. 4.8. Индексная матрица, соответствующая выбранным разбие-
нию и нумерациям, представлена таблицей.
Номер элемента Локальные номера
1 2 3
Глобальные номера
1 1 2 3
2 2 4 3
3 4 5 3
Рис. 4.8. Область, состоящая
из трех элементов
Соответствие между глобальными
и локальными обозначениями неиз-
вестных температур имеет вид:
ы,(,) = Hj, и^2) = «2 , н|3) = w3 ,
«1(2) =«2> «22) =м4» «32) =м3« 0-32)
,/3) _ „(3) -и
Теперь рассмотрим структуру гло-
бальной матрицы и глобального век-
тор-столбца. Начнем с первого урав-
76
нения. Поскольку узел 1 содержится только в первом элементе, то
в первом уравнении (4.21) остается только частная производная от
функционала первого элемента и оно принимает вид
az az(1) Q
du| duj
В соответствии с (4.27), (4.28) в локальных координатах это урав-
нение записывается в виде
„d^d) . „(Dud) . „(Dj-d) __ „(D
*11 “1 + «12u2 + *13U3 _<P1 •
Для первого элемента локальные и глобальные номера совпадают
(4.32), и поэтому окончательно первое уравнение системы записыва-
ется так:
„d)„ + „(•),. +„(Du — m(D
611 М1 + «12 “2 + «13 и3 - *П •
Отсюда вытекает, что являются первым, вторым
и третьим коэффициентами первой строки глобальной матрицы G,
a - первым коэффициентом глобального вектор-столбца:
G„=tS’. < Й’.О. о).
Сложнее обстоит дело со вторым узлом, принадлежащим двум
элементам 1 и 2. Второе уравнение системы имеет вид:
а^=а^=0
ди2 ди2
Из (4.27) и (4.28) вытекает, что в локальных обозначениях неиз-
вестных это уравнение записывается следующим образом:
-(Dyd) . _(1)U(D . „(Dud) _m(D +
«21 "1 ’522 “2 + «23и3 Ф2 +
+ </2М2) + „(2)и(2) + «/2М2) _ го(2) - О
+ 5Ц "1 +512 м2 + «!3 м3 Ф1
Обратим внимание, что при записи 5J(2)/du2 использованы коэф-
фициенты первой строки локальной матрицы и первый коэффициент
столбца для второго элемента, поскольку температура и2 в его ло-
кальных обозначениях имеет первый номер: u2 = Wi<2). Теперь для по-
лучения искомого уравнения необходимо в соответствии с (4.32) за-
менить локальные обозначения неизвестных на глобальные:
77
eWu +₽(1)w - m(1)4-
621 U1 + 622 u2 + 623 M3 *2 +
+ ₽(2)M +P(2)W +₽<2>и -ro<2>-0
M2+612 “4+6|3 M3 *Pj ~u-
(4.33)
Группируя слагаемые с одинаковыми неизвестными и складывая
постоянные коэффициенты, получаем вторую строку глобальной
матрицы G и второй коэффициент глобального вектор-столбца Ф
в виде
ф2 = .^ М2’
(4-34)
Аналогичным образом можно получить и остальные три урав-
нения.
Алгоритм формирования глобальных матрицы и вектор-столбца.
Полученные выражения (4.34) позволяют изложить принцип форми-
рования m-го уравнения глобальной системы. Это формирование це-
лесообразно проводить путем постепенного суммирования вкладов
от различных элементов. При машинной реализации перед началом
формирования массивы, в которых помещают глобальные G и Ф,
обнуляются, а затем к их текущим значениям постепенно добавляются
соответствующие коэффициенты локальных матриц и столбцов. Яс-
но, что вклад в т-е уравнение системы дадут только те элементы,
у которых в строке индексной матрицы имеется номер т. Если т-й
узел числится в локальной нумерации какого-либо из этих элементов
l-м (I - 1, 2 или 3), то будет использована l-я строка локальной мат-
рицы g(n)
и l-vi коэффициент локального вектор-столбца
. Найден-
ный нужный коэффициент локального столбца прибавляется к теку-
щему значению m-го коэффициента глобального столбца. Коэффици-
енты выделенной строки локальной матрицы элемента прибавляются к
соответствующим коэффициентам m-й строки глобальной матрицы,
имеющим порядковые номера, указанные в строке индексной матри-
цы, т.е. первому коэффициенту строки локальной матрицы соответ-
ствует первый номер «отсылки» в строке индексной матрицы, второ-
му коэффициенту - второй номер, третьему - третий.
Описанная процедура лежит в основе алгоритма формирования
глобальной матрицы и глобального вектор-столбца. Как было уже
отмечено выше, она реализуется путем последовательного перебора
элементов следующим образом. Берется первый элемент, анализиру-
ется его строка в индексной матрице и устанавливается, в какие три
уравнения этот элемент «дает вклад». Далее рассчитываются его ло-
кальная матрица и вектор-столбец. При этом расчете используется
78
информация о наличии у данного элемента граничных сторон, со-
держащихся в четвертом столбце индексной матрицы. Пусть локаль-
ным номерам 1, 2, 3 соответствуют фактические номера i, j, к. Тогда
первая строка локальной матрицы и первый коэффициент локально-
го вектор-столбца участвуют в формировании i-й строки глобальной
матрицы и i-го коэффициента глобального вектор-столбца. Произ-
водится сложение найденных локальных коэффициентов g^}\ g^,
g|^ с имеющимися значениями глобальных коэффициентов GiP G.„
(?rt. Затем аналогичная процедура повторяется для второй и третьей
строк локальной матрицы и второго и третьего коэффициентов ло-
кального столбца. Они участвуют в формировании строк глобальной
матрицы и коэффициентов глобального столбца с номерами j и к,
которые соответствуют локальным номерам 2 и 3.
Изложенный на примере треугольных элементов разбиения метод
формирования глобальных матрицы и вектор-столбца, основанный
на введении локальной нумерации узлов и неизвестных, легко пере-
носится и на случай более сложных элементов разбиения. Он является
наиболее общим, часто используемым и тем более эффективным, чем
сложнее применяемые конечные элементы.
Свойства системы разностных уравнений и методы ее решения.
Теперь рассмотрим ряд важных свойств, которыми обладает гло-
бальная матрица. Во-первых, можно доказать, что она является
симметричной. Во-вторых, глобальная матрица для задач большой
размерности М является сильно разреженной, т.е. большинство ее
элементов - нулевые. Наконец, путем введения разумной нумера-
ции узлов ее можно сделать ленточной.
Остановимся на двух последних свойствах матрицы. Очевидно,
что коэффициент Gmj в m-й строке глобальной матрицы отличен от
нуля, только когда узлы с номерами muj являются вершинами како-
го-то общего для них элемента. В этом случае в строке индексной
матрицы, соответствующей этому общему элементу, будут содержать-
ся номера т и j. Указанное обстоятельство и объясняет разреженность
глобальной матрицы, поскольку, например, для треугольных элемен-
тов при значительном числе треугольников N большинство возмож-
ных пар узлов т и к не являются вершинами общего треугольника,
и, следовательно, соответствующий элемент глобальной матрицы
Ч* = 0.
Рассмотрим влияние нумерации узлов на структуру глобальной
матрицы G. Из сказанного выше вытекает, что расположение нуле-
вых элементов в матрице зависит от способа нумерации узлов. На-
79
пример, в рассмотренном выше конкретном примере при нумерации,
указанной на рис. 4.8, глобальная матрица G выглядит так:
Сн
G2i
О
^12 G|3
G22 623
G32 G33
g42 g43
0 G53
0 0
G24 0
G34 g35
G44 G45
G54 Gi5
(435)
О
Если же перенумеровать узлы так, как это сделано на рис. 4.9, то
матрица G примет вид:
Gu
О
g12
g32
о
о
^23
^33
G43
^53
О <?15
О <?25
G34 <^35
G44 g45
G54 g55
(4.36)
_^51 G52
Общее число нулевых элементов в (4.36) не изменилось по сравне-
нию с (4.35), однако в (4.35) нулевые элементы расположены лишь на
главной диагонали и на двух прилегающих к ней верхних и нижних
диагоналях, а в (4.36) эти элементы «разбросаны» по всей матрице.
Таким образом, при разумной нумерации узлов глобальная мат-
рица G имеет ленточный вид, т.е. все
ненулевые коэффициенты расположе-
ны в пределах полосы, образованной
рядом верхних и нижних диагоналей,
5 9^ - примыкающих к главной диагонали.
/ X. / Из симметрии матрицы следует, что
/ у. (*' / число верхних и нижних диагоналей
/ (2) X. / с отличными от нуля коэффициентами
/ Х^/ одинаково.
♦. 2 Поскольку для треугольного раз-
биения коэффициент Gmk отличен от
Рис. 4.9. Переименование узлов НУЛЯ только в случае, когда узлы т и к
области из трех элементов принадлежат одному треугольнику, то
положение наиболее удаленного от
главной диагонали ненулевого элемента матрицы определяются мак-
симальной по всем парам общих вершин треугольников разностью
номеров узлов, т.е. величиной
80
L = max(n)
m —
(437)
где m,k- номера узлов л-го треугольника.
Схематичный вид глобальной ленточной матрицы показан на
рис. 4.10. Символом «х» обозначены ненулевые коэффициенты. Все
коэффициенты, расположенные за пределами полосы, ограниченной
линиями, параллельными главной диагонали, равны нулю. В общем
случае коэффициенты встречаются и внутри полосы. Число (L + 1)
будем называть шириной полосы (шириной ленты) матрицы.
Ленточный характер и симметрия глобальной матрицы позволя-
ют значительно сократить объем машинной памяти, требуемой для
ее хранения. При программировании задачи предусматривается за-
пись матрицы не в виде массива длиной М х М, а в виде массива, со-
держащего лишь элементы, находящиеся в пределах полосы на глав-
ной диагонали и выше.
Например, если требуется решить задачу с числом неизвестных
узловых температур М = 300, то для записи матрицы в общем виде
необходимо хранить 300 х 300 — 9 104 вещественных чисел. Пусть
ширина ленты в этой задаче равна L + 1 = 40. Тогда запись матрицы
в сокращенном виде потребует массив длиной (L+1)M - 12 104 элемен-
тов, если матрицу запом-
нить в виде «прямоуголь-
ника» с £ + 1 столбцами
и М строками, т.е. для об-
легчения логики програм-
мирования предусматри-
вать место для хранения
фиктивных элементов в пос-
ледних L строках (зашт-
рихованный треугольник
вне матрицы на рис. 4.10).
Если же не учитывать фик-
тивные элементы «хвоста»,
то потребуется запомнить
[(£+1>И -£(£+1)/2] =
= 1,12 -10 4 чисел.
Рис. 4.10. Схематичный вид глобальной
ленточной матрицы
Наконец перейдем к вопросу решения системы уравнений. Для
решения систем уравнений МКЭ применяются как прямые, так и ите-
рационные методы. Причем последние обычно используют в тех слу-
чаях, когда объем оперативной памяти не позволяет хранить всю
81
глобальную матрицу даже с учетом ленточного симметричного вида.
Из прямых методов хорошо зарекомендовал себя на практике и по-
лучил широкое распространение метод квадратного корня. Этот ме-
тод пригоден только для систем линейных уравнений с симметрич-
ной матрицей и по затратам машинного времени примерно вдвое
быстрее метода исключения Гаусса. В математическом обеспечении
ЭВМ имеются стандартные программы, реализующие метод квад-
ратного корня. Предусмотрен и случай систем с ленточной матрицей.
В заключение подчеркнем, что использование той или иной стан-
дартной подпрограммы решения системы уравнений требует опреде-
ленного способа записи глобальной матрицы в одномерный массив.
Применяемые способы различны для разных подпрограмм, т.е. мо-
жет ограничиваться запись по строкам или по столбцам, с учетом
или без учета «хвоста» ленты. Это обстоятельство следует учитывать
при программировании алгоритма формирования глобальной мат-
рицы.
4.5. Программная реализация МКЭ
Существенным достоинством МКЭ является возможность состав-
ления программ численного расчета полей в областях сложной гео-
метрической конфигурации, которые проще по логической структуре
и по заданию исходных данных, чем программы, реализующие метод
конечных разностей для таких областей. В данном подразделе рас-
смотрим в качестве примера структуру программы для решения дву-
мерной задачи (4.1), (4.2) в областях произвольной формы при тре-
угольных элементах разбиения.
В программе решения задачи методом конечных элементов вы-
полняются следующие основные процедуры.
1. Разбиение области на элементы, нумерация элементов, гло-
бальная и локальная нумерации узлов и формирование на их основе
индексной матрицы.
2. Формирование глобальных матрицы и вектор-столбца системы
алгебраических уравнений, реализуемое на основе расчета локаль-
ных матриц и столбцов.
3. Решение системы разностных уравнений.
4. Расчет температур и тепловых потоков в различных точках
элементов разбиения, проводимый на основе принятой аппроксима-
ции температурного поля в элементе.
82
Остановимся на особенностях программной реализации первых
двух процедур, рассматривая их применительно к решению задач
(4.1), (4.2).
Автоматизация разбиения области. Простейший (но наиболее тру-
доемкий) способ реализации первой процедуры состоит в ручном
разбиении области D на треугольные элементы, ручной нумерации
узлов и дальнейшем вводе в качестве исходных данных массивов ко-
ординат узлов {хт , {ут и индексной матрицы. Однако в ре-
альных двумерных (и тем более трехмерных) задачах число узлов
и элементов может составлять несколько сотен, а иногда и тысяч,
и поэтому построение расчетной сетки вручную и ввод больших мас-
сивов чисел в качестве исходных данных нецелесообразны из-за зна-
чительных затрат времени на их подготовку и большой вероятности
появления ошибок. Следовательно, возникает задача автоматизации
процедуры разбиения областей на элементы, нумерации элементов
и узлов и формирования индексной матрицы. При этом требуется
в качестве входной ин
рмации для соответствующей подпрограм-
мы задавать сравнительно небольшое число данных, описывающих
геометрию области сложной формы и густоту сетки, а на ее выходе
получать массивы координат узлов и индексной матрицы, позво-
ляющие перейти к формированию глобальной матрицы системы
уравнений МКЭ.
Для областей типа прямоугольника, треугольника или некоторых
более сложных, но конкретных конфигураций можно составить ча-
стные подпрограммы, реализующие заполнение массивов координат
узлов и индексной матрицы на основе данных о размерах области
и шагах по осям х, у. Например, можно покрыть область прямо-
угольной сеткой, а затем построить треугольники путем разбиения
элементарных прямоугольников диагоналями, как это сделано на
рис. 4.11. Программирование таких процедур сводится к организа-
ции нескольких циклов с переменными пределами.
Однако для практических приложений большой интерес пред-
ставляют универсальные программы автоматического разбиения об-
ластей различной сложной формы. В литературе предложены разные
способы описания геометрии и алгоритмы дискретизации областей
для МКЭ. Наибольшее распространение получил следующий подход.
Область сложной формы разбивается вручную на подобласти, кото-
рые называются «макроэлементами».
Эти подобласти должны достаточно хорошо описывать геомет-
рию расчетной области. Обычно макроэлементы выбирают в форме
треугольников и выпуклых четырехугольников, но иногда использу-
83
ют и подобласти, ограниченные кривыми второго порядка. Число
таких макроэлементов обычно невелико (несколько единиц или де-
сятков для сложных областей), и поэтому это разбиение можно опи-
сать путем задания координат узлов макроэлементов и некоторой
условной нумерации макроэлементов и их узлов.
Затем реализуется автоматическое разбиение каждого из макро-
элементов на элементарные треугольные элементы. Для этого в ис-
ходных данных лишь указывают параметры, характеризующие гус-
тоту сетки в каждом из макроэлементов.
Приведем два примера автоматизации разбиения. В первом ис-
пользуются макроэлементы в виде выпуклых четырехугольников
(рис. 4.12). Для каждого макрочетырехугольника задается следую-
щая информация: номера узлов в вершинах и признаки принадлеж-
ности сторон границе, числа дроблений на отрезки по двум смежным
ии
сторонам, массивы коор,
нгат вершин четырехугольников. Каждый
из макрочетырехугольников автоматически разбивается на малень-
кие четырехугольники прямыми, проходящими через точки разбие-
ния противоположных сторон на равные отрезки.
1 2 4 7 х
Рис. 4.11. Разбиение прямоугольных
элементов на треугольные
Рис. 4.12. Представление области
разбиения макроэлементами в виде
выпуклых четырехугольников
Полученные четырехугольники, в свою очередь, разбиваются ко-
роткими диагоналями на элементарные треугольники.
В результате число треугольников в каждом макроэлементе равно
2 к\к2 , где к\, к2 ~ числа отрезков на смежных сторонах. Заметим, что
при таком способе для общих сторон макроэлементов должна зада-
ваться одинаковая кратность дробления.
84
Другой способ автоматизации разбиения иллюстрирует рис. 4.13.
Здесь в качестве макроэлементов взяты треугольники. Информация
о них задается почти в таком же виде, как и для «элементарных» тре-
угольников: массивы координат вершин и индексная матрица, но
с одним отличием. В строке индексной матрицы для каждого макро-
элемента содержится еще одно число - кратность дробления к. Если
к = 0, то макроэлемент не дробится и принимается в качестве конеч-
ного элемента. При к = 1 путем соединения центров сторон прово-
дится разбиение макроэлемента на четыре подобных треугольника
(рис. 4.13). При к = 2 каждый из полученных четырех треугольников
еще раз разбивается на четыре подобных и т.д., т.е. число получен-
ных из макроэлемента треугольников равно 4*. Кратность дробления
соседних макроэлементов может различаться не более чем на едини-
цу. При этом, чтобы избежать появления «лишних» узлов на стороне
треугольника с меньшей кратностью дробления, автоматически прово-
дится построение еще нескольких треугольников. Для этого узел, ле-
жащий на стороне треугольника, соединяется с противоположной
вершиной, как это показано на рис. 4.13 пунктирными линиями. Дос-
тоинством данного способа разбиения является возможность резко
сгущать сетку в областях с большими градиентами температур, ис-
пользуя при этом сравнительно небольшое число макроэлементов.
Приведенные на «словесном» уровне описания методик автома-
тического разбиения двумерной области могут создать впечатление
о простоте их программной реализации. Однако на самом деле про-
граммное «воплощение» этих алгоритмов требует довольно значи-
тельных усилий. к -2
В процессе автоматического раз- .Л.
биения области проводится и нумера- / X/ Х/ У у*
ция получившихся элементов разбие- . _
ния. После этого возникает задача К /\/7
оптимальной перенумерации узлов X/
с целью уменьшения ширины ленты
матрицы. В областях, имеющих форму, \
близкую к прямоугольной, нумерацию к = 1 \ / \ д
узлов целесообразно проводить по- \\
следовательно вдоль направлений, ^******^^
параллельных короткой стороне
и содержащих меньшее число узлов. ?ис- 4-13- Разбиение треуголь-
ников на более мелкие
85
12
Рис. 4.14. Нумерация узлов тре-
угольников области разбиения
В некоторых программах нумерацию проводят последовательно
в пределах каждого макроэлемента в порядке их обхода. При этом
в треугольниках, лежащих у границ макроэлементов, могут возни-
кать большие разности номеров узлов. Можно предложить и другие
подходы. Например, хорошо зарекомендовал себя для решения мно-
гих практических задач следующий эвристический способ перенуме-
рации узлов. Узел с номером 1 выбирают где-либо на границе облас-
ти. Номера 2,3 и т.д. присваивают узлам, которые соединены с узлом
номер 1 стороной элементарного треугольника (рис. 4.14). После то-
го как пронумерованы все узлы,
соединяющиеся с узлом номер 7,
в качестве «опорного» берется
узел номер 2 и продолжается ну-
мерация узлов, соединенных
с ним и не пронумерованных ра-
нее. Подобная процедура последо-
вательной нумерации узлов, соеди-
ненных с наименьшим из ранее
присвоенных номеров, продолжа-
ется до полной нумерации всех
узлов сетки. После такой перену-
мерации подсчитывается полу-
чившаяся ширина ленты матрицы,
которая также является выходным параметром программы разбие-
ния и перенумерации. Поскольку программная реализация описан-
ных автоматических процедур довольно сложна, примеры соответст-
вующих подпрограмм здесь не приводятся.
Формирование глобальных матрицы и вектор-столбца. Перейдем
к рассмотрению реализации алгоритмов формирования глобальных
матрицы и вектор-столбцов. Обычно они строятся на основе методи-
ки, изложенной в предыдущем параграфе. В этом случае оформляется
подпрограмма для расчета локальных матриц и столбцов.
Входными данными для нее являются координаты узловых точек,
передаваемые в соответствии с локальной нумерацией узлов, инфор-
мация о сторонах элемента, лежащих на границе области, и значения
параметров Хх, X, qv, qs, а . Последние могут вычисляться с помощью
специальных подпрограмм, в которых входными параметрами явля-
ются координаты узлов элемента, а выходными - значения Хх, X, qv,
qt, а для этого элемента. В приводимом ниже примере предусмотрено
использование двух таких подпрограмм: для расчета Хх, Х^, q»
(подпрограмма ААА) и qs, а (подпрограмма ВВВ).
86
1 SUBROUTINE MKE (X, Y, IND, M, N, AAA, BBB, MS, G, F)
2 C
3 С ПОДПРОГРАММА ФОРМИРОВАНИЯ МАТРИЦЫ МЕТОДА
КОНЕЧНЫХ
4 С ЭЛЕМЕНТОВ ДЛЯ ДВУХМЕРНОЙ, СТАЦИОНАРНОЙ ЗАДАЧИ
ТЕПЛОПР.
5 С
6 С ОПИСАНИЕ ПАРАМЕТРОВ:
7 С М-ЧИСЛО УЗЛОВ
8 С
9 С
10 С
11 С
12 С
13 С
14 С
15 С
N - ЧИСЛО КОНЕЧНЫХ ЭЛЕМЕНТОВ
Х(М), Y(M) - МАССИВЫ КООРДИНАТ УЗЛОВ
IND (4*N) - ИНДЕКСНАЯ МАТРИЦА
AAA, ВВВ - ИМЕНА ВНЕШНИХ ПОДПРОГРАММ
MS - ШИРИНА ЛЕНТЫ ЛЕНТОЧНОЙ МАТРИЦЫ
G - ВЫХОДНАЯ ЛЕНТОЧНАЯ МАТРИЦА
F - ВЫХОДНОЙ ВЕКТОР-СТОЛБЕЦ ПРАВЫХ ЧАСТЕЙ
ВХОДНЫЕ
ПАРАМЕТРЫ
ВЫХОДНЫЕ
ПАРАМЕТРЫ
16 DIMENSION X(l), Y(l), IND(l), G(l), F( 1)
17 С
18 С ВЫЧИСЛЕНИЕ ШИРИНЫ ЛЕНТЫ MS
19 С
20 MS = 0
21 LL —0
22 DO 1 L = 1,N
23 I = IND (LL)
24 I = IND (LL + 1)
25 К = IND (LL + 2)
26 MS = MAXO (MS, IABS (I - I)
27 *IABS (I - K), IABS (I - K)
28 1 LL = LL + 4
29 C
30 С ОЧИСТКА МАССИВОВ G, F
31 C
32 I = M + MS • (2*M - MS - 1) / 2
33 DO 2 I = 1,1
34 2G(l) = 0
35 DO3 1= 1,M
36 3 F (I) = 0
37 C
38 С НАЧАЛО ЦИКЛА ПО ТРЕУГОЛЬНИКАМ
39 С
40 LL = I
41 DO 12 L- 1,N
42 C
43 С НОМЕРА УЗЛОВ
44 С
45 I = IND (LL)
46 I = IND (LL + 1)
47 К = IND ( LL + 2)
48 C
49 С ВЫЧИСЛЕНИЕ КООРДИНАТ ЦЕНТРА ТРЕУГОЛЬНИКА
50 С
87
51
52
53 С
54 С
55 С
56
57
58
59
60
61
62 С
63 С
64С
65
66
67
68
69
70
71 С
72 С
73 С
74
75
76 С
77 С
78 С
79 С
80
81
82
83
84
85
86
87
88
89
90
91
92 С
93 С
94 С
95 С
96
97
98
99 С
100 С
101 С
102
ХС = (Х(1) + Х(1) + Х(К)) / 3
YC = Y(I) + Y(I) + Y(K)) / 3
ПРИВЕДЕНИЕ КООРДИНАТ УЗЛОВ К ЦЕНТРУ
XI - Х(1) - ХС
Х2 = Х(1) - ХС
ХЗ = Х(К) - ХС
Yl = Y(I) - YC
Y2 = Y(I) - YC
Y3 = Y(K) - YC
ВЫЧИСЛЕНИЕ КОЭ
ИЦИЕНТОВ ФУНКЦИИ ФОРМЫ
Bl — Y2 - Y3
С1 =ХЗ-Х2
В2 = Y3-Y1
С2 = XI - ХЗ
ВЗ = Yl - Y2
СЗ = Х2-Х1
ВЫЧИСЛЕНИЕ ПЛОЩАДИ ТРЕУГОЛЬНИКА
А = (X2*Y3 - X3*Y2 + XI *Y2 -
X1*Y3 + X3*Y1 - X2*Y1) / 2
ВЫЧИСЛЕНИЕ ЭЛЕМЕНТОВ ЛОКАЛЬНОЙ МАТРИЦЫ
И ВЕКТОРА ПРАВЫХ ЧАСТЕЙ
CALL AAA (ALX, ALY» ALFV, QV, XC, YC)
H = 4* A
D = ALFV*A / 6
Gil = (ALX*B1*B1 + ALY*C1*C1)/H + D
G12 = (ALX*B1*B2 + ALY*C1*C2)/ H + D / 2
G13 = (ALX*B1*B3 + ALY*C1*C3) / H + D / 2
G22 = (ALX*B2*B2 + ALY*C2*C2) / H + D
G23 = (ALX*B2*B3 + ALY*C2»C3) / H + D / 2
G33 = (ALX*B3*B3 + ALY*C3*C3) / H + D
Fl = QV»A / 3
F2 = Fl
F3 = FI
ВЫЧИСЛЕНИЕ ДОПОЛНИТЕЛЬНЫХ СЛАГАЕМЫХ
ДЛЯ ГРАНИЧНЫХ ЭЛЕМЕНТОВ
LL = LL + 3
IF (IND (LL).EQ.O) GO ТО 5
IF (IND (LL).EQ.0) GO TO 4
ВКЛАД СТОРОНЫ IK
CALL BBB (ALFL. QL, (X(I) + X(K)) / 2, (Y(I) + Y(K)) / 2)
88
103
104
105
106
107
108
109
110
111 С
112 С
113 С
114
115
116
117
118
119
120
121
122
123 С
124 С
125 С
126 С
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153 С
154 С
Н = SQRT ((Х2 - Х3)**2 + (Y2 - У3)*ф2)
D- ALFL*H/3
С22 = С22 + D
СЗЗ = СЗЗ + D
C23 = C23 + D/2
D - QL ♦ Н / 2
Fl = Fl + D
F3 = F3 + D
ВКЛАД СТОРОНЫ IJ
4 CALL BBBfALFL, QL, (X(I) + X(I)) /2, (Y(I) + Y(I)) / 2)
H = SQRT ((XI - X2)**2 + (Y1 - Y2)**2)
D = ALFL*H / 3
Gil =C11 + D
G22 = G22 + D
G12 = G12 + D/2
D = QL * H / 2
Fl = Fl +D
F2 = F2 + D
ДОБАВЛЕНИЕ ЭЛЕМЕНТОВ ЛОКАЛЬНОЙ
МАТРИЦЫ К ЭЛЕМЕНТАМ ГЛОБАЛЬНОЙ
5N1 =1
N2 = I
D-G11
ASSIGN 6 ТО LBL
GOTO 13
6 N2 = I
D = G12 .
ASSIGN 7 TO LBL
GOTO 13
7N2 = K
D = GI3
ASSIGN 8 TO LBL
GO TO 13
8N1 = 1
N2 = I
D = G22
ASSIGN 9 TO LBL
GOTO 13
9 N2 = К
D = G23
ASSIGN 10 TO LBL
GOTO 13
10N1 =K
D = G33
ASSIGN 11 TO LBL
GOTO 13
TO ЖЕ ДЛЯ ВЕКТОРА ПРАВЫХ ЧАСТЕЙ
89
155 С
156 11 F(I) = F(I) + Fl
157 F(I) - F(I) + F2
158 F(K) = F(K) + F3
159 C
160 С КОНЕЦ ЦИКЛА ПО ТРЕУГОЛЬНИКАМ
161 С
162 12 LL = LL 4-1
163 RETURN
164 С
165 С ФОРМИРОВАНИЕ МАТРИЦЫ В ВИДЕ, ПРИГОДНОМ
166 С ДЛЯ ОБРАЩЕНИЯ К ПОДПРОГРАММЕ МСНВ
167 С
168 13 Ml = MINO (Nl, N2)
169 М2 - МАХО (Nl, N2)
170 ММ=1+М2-М1
171 MS1 = MS+1
172 MMS = M-MS
173 IF (Ml, GE, MMS'i GOTO 14
174 NN = MS1*(M1-1) + MM
175 GOTO 15
176 14 МММ = M-Kfl + 1
177 NN = MS1*MMS + MS*MS1 / 2 - MMM*(MMM+1) / 2 + MM
178 15 G(NN) = G(NN) + D
179 GOTO LBL (6,7, 8, 9,10,11)
180 END
Обращения к этим подпрограммам имеют вид:
CALL AAA (ХС, YC, АХ, AY, QV )
и CALL ВВВ (ХС, YC, ALF, QS).
Эти подпрограммы должны составляться пользователем для кон-
кретной задачи. При формировании матрицы происходит обращение
к ним с текущими координатами ХС, YC центра элемента (для ААА)
или центра стороны элемента (для ВВВ). Напомним, что при выводе
уравнений МКЭ мы считали свойства и мощности постоянными
в пределах элемента, поэтому в случае разрывных функций жела-
тельно, чтобы линии разрыва совпадали с границами элементов.
Входными данными для подпрограммы формирования глобаль-
ных матрицы и столбца являются: N - число треугольных элементов,
М - число узлов, MS - ширина ленты матрицы, X, Y-массивы коор-
динат узлов хт, ут длиной М, IND - индексная матрица - массив
длиной 4*N. Все эти данные получаются в результате выполнения
процедуры разбиения и перенумерации узлов.
В индексной матрице для каждого треугольного элемента (и = 1, ...,
N) записываются четыре целых числа. Первые три числа соответству-
ют глобальным номерам /, j, к узлов данного треугольника, записан-
ным в порядке обхода против часовой стрелки, начиная с любой
90
вершины. Четвертое число - это признак, с помощью которого зада-
ется внешняя граница области. Напомним, что согласно принятому
нами способу этот признак принимает значения: 0 - если ни одна из
сторон треугольника не принадлежит границе, 1 - если граничной
являются две стороны между узлами i и j, 2- если граничной являют-
ся две стороны между узлами i и j и между узлами j и к. Другие си-
туации с расположением граничных сторон можно исключить соот-
ветствующим выбором порядка записи номеров i,j, к.
В приводимой подпрограмме из-за малой размерности локальных
матрицы и столбца для расчета их коэффициентов не используется
специальная подпрограмма, а они вычисляются непосредственно
в ней по формулам (4.29), (4.30).
В подпрограмме производятся следующие процедуры. Сначала
обнуляются массивы G и F, соответствующие глобальным матрице
и столбцу. Далее организуется цикл по всем элементам, в котором:
• по индексной матрице определяются глобальные номера узлов /,
/, к в данном элементе и принадлежность сторон границе;
• вычисляются коэффициенты функций формы элемента по (4.14)
и рассчитываются коэффициенты локальных матрицы и столбца по
(4.29), (4.30);
• значения всех коэффициентов локальных матрицы и вектор-
столбца прибавляются к текущим значениям коэффициентов гло-
бальных матрицы и вектор-столбца, индексы которых определяются
глобальными номерами узлов i,j, к .
В результате обхода всех N треугольников оказываются полно-
стью сформированными глобальная матрица G и столбец правых
частей Ф системы уравнений.
Отметим, что при формировании матрицы G необходимо учиты-
вать способ записи матрицы в машинной памяти для используемой
стандартной подпрограммы решения системы линейных уравнений.
При этом коэффициенты матрицы должны быть записаны в одномер-
ный массив путем последовательного обхода верхней части ленты над
главной диагональю по строкам. Такой пересчет индексов элемента мат-
рицы в индекс одномерного массива реализован операторами 168-177.
Контрольные вопросы
1. В чем суть метода конечных элементов для решения задач теплопро-
водности?
2. Как производится построение дискретной модели и функции формы
элементов?
3. Как составляется система уравнений метода конечных элементов?
91
Глава 5
МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ
ДЛЯ РЕШЕНИЯ ЗАДАЧ
ПЛАСТИЧЕСКОГО
ФОРМОИЗМЕНЕНИЯ
Метод конечных элементов (МКЭ) включает различные подходы,
в которых для определения напряжения, деформации и перемещения
материал условно разбивают на отдельные элементы, соединенные
в узловых точках. Этот метод впервые был использован при упругом
анализе конструкций, а затем нашел применение для решения самых
разнообразных инженерных проблем, в том числе для исследования
процессов пластического формоизменения.
В данной главе, следуя [4], рассмотрено применение МКЭ к упру-
гопластическим и жесткопластическим материалам, а также сумми-
рованы результаты анализа процессов формоизменения с использо-
ванием МКЭ.
5.1. Основные характеристики элементов
5.1.1. Конечные элементы
Представим, что заготовка в состоянии плоской деформации раз-
делена на конечное число треугольных призматических элементов
с воображаемыми границами (рис. 5.1); элементы соединены только в
узловых точках (узлах) и силы не могут передаваться через боковую
поверхность элементов. Элементы и узлы так пронумерованы, что
соседние элементы и узлы имеют близкие номера.
При небольшом перемещении инструмента узлы перемещаются
в новые положения. Перемещения узлов считаются неизвестными па-
раметрами, которые позднее предстоит определить при заданных
граничных условиях. В двумерной модели неизвестных параметров
вдвое больше, чем узлов, поскольку каждому узлу соответствуют две
компоненты перемещения и со> в направлениях х и у. Приращения
напряжении и деформаций могут быть определены в результате вы-
числения перемещения узлов. Нелинейность, связанная с геометриче-
92
13
18
14
Рис. 5.1. Конечные элементы и индек-
сация (нумерация) узлов и элементов
скими изменениями, а также со свойствами материала при пластиче-
ской деформации, для проведения вычислений требует разбить весь
процесс деформирования на ряд этапов. Для получения соответст-
вующих значений координат, напряжений и деформаций вычислен-
ные приращения следует прибавлять к значениям, полученным на
предыдущих этапах.
Изложенный выше метод,
основанный на рассмотрении
смещения узлов, называется ме-
тодом перемещений. Существует
метод, в котором основными
неизвестными параметрами яв-
ляются напряжения: это метод
сил (напряжений). Читатель
может догадаться, что первый
метод имеет много общего с мето-
дом верхней оценки, при котором
рассматривают кинематически
возможные поля перемещений
(или скоростей), а второй - с ме-
тодом нижней оценки, в котором
используют статически возмож-
ные поля напряжений. Основ-
ным недостатком метода сил
является необходимость учета
большого числа ограничивающих условий. С другой стороны, боль-
шая степень свободы, присущая методу смещений, позволяет полу-
чить различные результаты при решении задач упругопластической
деформации. '
10
5.1.2. Функции перемещения
В качестве примера рассмотрим треугольный элемент ijm в усло-
виях плоской деформации (рис. 5.2). Толщина элемента S, координа-
ты Эйлера-Декарта трех узлов (узловых точек) соответственно (х, у.),
(х, ), (хт, уп). При небольшом изменении шага деформации узлы
перемещаются в направлениях х и у за время Az на (сох., со ), (со , со ),
(солт, со Если смещение узла невелико и дальнейшая деформация
происходит непрерывно, то компоненты скоростей узлов (и., сх), (и, сх),
(ит, uj можно определить делением перемещения на Az. Для удобства
93
рассмотрения вместо перемещений узлов обычно используют их ско-
рости.
Компоненты перемещения или скорости элемента интерполиру-
ются функцией перемещения. Простейшими функциями перемещения
являются линейные полиномы
u = Oj + d2x + d3x;
и = а4 + а5х + а6у.
(5.1)
а1; а2,...,а6 определяются путем подстановки в эти уравнения ско-
ростей и координат узловых точек, т.е.
Рис. 5.2. Треугольный элемент
с тремя узлами
ut = cq + а2х, + а3у;
Uj = dj + a.2Xj + а3уу; (5.2)
= а1 + <X2*m + «3?т-
Решение этих уравнений по-
зволяет найти
где А - площадь треугольника
ijm, т.е. А = (ххут + хту{ + х,уу -
Х}У< хгУт хтУ j
а а, Ь., с определены следующими выражениями:
= ^]Ут ~~ ХтУjr
~~ У jm ~~ У / ” Ут ’
г. = y . = X — X •
mj т ) ’
Другие коэффициенты получаются путем циклической переста-
новки индексов. Используя эти значения, уравнение (5.1) можно пе-
реписать в терминах скоростей узлов:
« = { (а< + btx + с, у + (aj +Ь}х + с/У)и j+(am+ bmx + сту )ит }/2 А;
и = {(а, + Ь(х + c.yjo,- + (а7 + bj-x + Cjy})j + (am + bmx + cmy\>m]/2A.
94
Если используют функцию линейного перемещения согласно
уравнению (5.1), то требование непрерывности скорости через гра-
ницу, проходящую между двумя смежными элементами, удовлетво-
ряется автоматически, поскольку смещение изменяется непрерывно
вдоль границы между двумя узловыми точками в обе стороны.
5.1.3. Скорость деформации и вращение
В задачах плоской деформации компоненты скорости малых де-
формаций определяются как
ё v = ди / дх;
•Л
ёу=ди/ду; /
уху = ди/ду+ дх.
При подстановке уравнения (5.4) в (5.5) компоненты скорости де-
формации треугольного элемента выражаются через скорости узлов
Вообще говоря, эти уравнения можно переписать в матричной
форме:
{ё} = №}, (5.7)
где [ ] и {} обозначают матрицы, содержащие по единственному столбцу,
т.е. вектор. Матрица [2?] изменяется в зависимости от типа элемента и
функции перемещения. При плоской деформации треугольного эле-
мента точные выражения для {ё},р J {«} следующие:
Г«. 1
о
bj
О
CJ
О
CJ
bj
Ьщ
Поскольку скорости деформации не содержат слагаемых с коор-
динатами х и у, скорость деформации постоянна по всему элементу.
95
Скорость деформации изменения объема su при плоской дефор-
мации
£о Ёх (5.9)
Для треугольного элемента
£и — iPiUj + Cj-D( + bjllj + Cjpj + ^nt^m ‘ (5. Ю)
Скорость изменения объема элемента (объем V = AS)
V = kuV = ev(AS).
(5.П)
Угловая скорость вращения элемента против часовой стрелки
coXJ, = (ди/дх-ди1ду)/2.
(5.12)
При плоской деформации треугольного элемента угловая ско-
рость вращения выражается через скорости узлов:
~ CjUj ~~CmUm (5.13)
5.1.4. Степень свободы деформации
Степень свободы деформации элемента ограничена условием не-
сжимаемости при пластической деформации. Для элемента с малой
степенью свободы действительную деформацию нельзя представить
соответствующим образом, и вычисление нагрузки приводит обычно
к завышенным результатам.
Степень свободы элемента в системе с большим числом элементов
определяется как произведение степени свободы узла (две при пло-
ской деформации) на число узлов в элементе. В большой области,
заполненной треугольными элементами, число элементов вдвое
больше числа узлов, т.е. число узлов, приходящихся на один элемент,
равно 1/2, таким образом, степень свободы элемента N = 1. Когда
накладывается условие несжимаемости, как в случае пластического
деформирования обычных металлов, степень свободы уменьшается.
Степень ограничения условием несжимаемости определяется числом
уравнений, выражающих условие несжимаемости элемента.
С учетом уравнения (5.10) условие несжимаемости треугольного
элемента можно записать в виде
— ^"m^?n)/24 — 0. (5.14)
96
Если это уравнение выполняется в точке, то оно удовлетворяется
всюду в элементе, таким образом, число ограничений Nc — 1.
Для определения степени свободы при условии несжимаемости
в качестве критерия принимают отношение Nj /Nc . Когда это отно-
шение намного больше единицы, считается, что деформация заклю-
чена в элементе. Если отношение меньше единицы, необходимо смяг-
чить условие несжимаемости. Для треугольного элемента в условиях
пластической деформации NfjNc=\. Эмпирически при анализе про-
стой деформации, такой как в условиях испытания на растяжение,
треугольный элемент проблем не создает, но для рассмотрения слож-
ной деформации при формоизменении металла не всегда удобен.
5.1.5. Характеристики некоторых элементов
А. Элемент плоской деформации. На рис. 5.3 показаны два четы-
рехугольных элемента с различными функциями перемещения.
Рис. 5.3. Четырехугольные элементы с четырьмя узлами:
а - изо пара метрический четырехугольник; 6 - пересеченный
диагональю четырехугольник
Поскольку четырехугольный элемент имеет четыре узла, скорость
нельзя представить одной единственной линейной функцией узла.
Для интерполяции скорости необходимо использовать некоторые
функции перемещения.
В случае изопараметрического четырехугольного элемента
(рис. 5.3, а) координаты Эйлера-Декарта переходят в новую коорди-
натную систему £ и т). Новые координаты выбраны так, что линии
постоянных Е, и и являются прямыми в плоскостях х - у и принимают
97
значения ±1 на сторонах четырехугольника. Интерполяционные
формулы двух систем координат представляются двумя уравнениями
с координатами четырех узлов (х., у), (xjf у.), (хж, yj и (хк, ук).
х={(1 - - nk + О+ЭД - п)*,-+(i+ЭД+пК, + (1 - ЭД+пК }/4,
У = {(1 - ЭД - п)л + (1 + ЭД - п).Уу + 0 + ЭД+ + (1 - ЭД+п)№ }/4 •
Те же формулы можно использовать для определения компонент
скорости элемента:
ч={(1 - ЭД> -nk- + (1+ЭД1 - nk7 + О+ЭД1+т&т + 0 - ЭД+ЭД* }/4;
v = j(i - ЭД - ЭД, +(1+ЭД-ЭД, ЭД+ЭД+ЭД», +(1-ЭД+ЭД*)/4-
Поскольку параметры E, и rj обычно используют для обозначения
координат и компонент скорости, элемент называют изопараметри-
ческим.
Компоненты скорости деформации записывают в форме уравне-
ния (5.7) как линейные функции восьми компонент скоростей (и., и
и„ о „ и , и , и., и.).
У’ У’ т т к7 к'
Матрицы [В] есть
д
о
о
С,
Д
Bj о
О Cj
cj BJ
вт 0 вк
О ст о
(5.17)
где J = {{-yim + ук.) + (у. + у^ {(-х* + хтк) +(х.. + хтк) п) - {(-*>, +
+ХШ > + <Х.7 + ^~Уц +У^ + +
-ХЛТ|
98
Степень свободы элемента N= 2, тогда условие несжимаемости
можно задать в виде
Ёи — (21 — yfofUj + XfoVj — у^ит + х^мт +
+ Уйп^к ~ ^йп^к ) Утк^( + ^mk^i Утк^] ~ ^mk^J +
+ УдМт ~ Xjj\)m — + XjjОд. уj^Uj + X^Vj +
4" У/k^j Xik^j ~ УцсЦщ + х^хлт + у— Х^\)^ )n}
= Di + D& + £>3T] = 0. (5.18)
Три уравнения = D2 = D3 = 0, чтобы обеспечить несжимаемость
по всей площади элемента, таким образом, Nc = 3 и N, /Nc = 2/3.
Четырехугольный элемент, представленный на рис. 5.3, б, состоит
из четырех треугольных подэлементов с внутренним узлом в точке
пересечения диагоналей. В каждом подэлементе поле скоростей ли-
нейно. Четырехугольный элемент обладает особым свойством при
плоской деформации, а именно, когда три из треугольных подэле-
ментов сохраняют свой объем постоянным, объем оставшегося, чет-
вертого, автоматически постоянен, таким образом, Nc = 3. Посколь-
ку степень свободы деформирования элемента равна 4 (2 для четы-
рехугольника и 2 для центрального узла), TVy /Nc -4/3.
Значения Nf,Nc и Nj/Nc приведены в табл. 5.1 для элементов,
испытывающих плоскую д
рмацию, а также для осесимметричных
элементов.
Значения степеней свободы для различных элементов
Таблица 5.1
Тип элемента Деформация
плоская осесимметричная
"г N( N<
Треугольный Изопараметрический 1 1 1 1 3 1/3
четырехсторонний 2 3 2/3 2 5 2/5
Поперечный четырех- сторонний 4 3 4/3 4 12 1/3
Для плоской деформации Nf /Nc максимально для поперечного
четырехугольника и минимально для изопараметрического четырех-
угольника. Для точного анализа задач плоской деформации жела-
тельно использовать поперечный четырехугольный элемент. При
решении практических задач изопараметрический элемент также
99
можно использовать, если ограничение смягчить оценкой скорости
деформации изменения объема в одной или в небольшом числе гаус-
совых точек элемента.
Б. Элементы для плоского напряженного состояния. Все элементы,
применяемые для анализа плоской деформации, можно использовать
и в задачах плоского напряженного состояния. Свободное деформи-
рование в направлении толщины смягчает требование несжимаемо-
сти, и выбор элемента не оказывает серьезного влияния на схему де-
формации.
В. Осесимметричные элементы. Для анализа осесимметричного
напряженного состояния используют элементы, имеющие форму
кольца с треугольным или четырехугольным меридиональным сече-
нием (рис. 5.4).
Функции перемещения в точности те же, что и для элементов при
плоской деформации, только х, у, z заменяют соответственно на
эйлеровы цилиндрические координаты г, z, <р.
Компоненты скорости малой деформации:
ёг-ди1дг\ Ё9=и1г;
e^dvldz; у_-ди! dz + dvl дг.
Л» * 9 9 &
(5-19)
Для треугольного кольцевого элемента
О bj 0 Ьт
О dj 0 dm
Cj 0 Cj О
Cj
— z
i Im
0
где a. = r. z -r z,,
l J m m f
+ b. + c z/r ,...
4
(5.20)
Компоненты скорости деформации ёг, ё2, yrz постоянны в эле-
менте, в то время как ёф содержит члены 1/r, z/r. Чтобы удовлетво-
рить условию несжимаемости, всюду в элементе необходимо выпол-
нить требование = Z>2 = £>з = 0, вытекающее из условия
ёг+ёф+ё2 = D{ + D'2 / г + D'3z / г -0,
таким образом, Nf!Nc = 1/3 (табл. 5.1).
100
г
6
Рис. 5.4. Осесимметричные элементы:
а - кольцо с треугольным сечением; б ~ кольцо с четырехугольным сечением
а
Поскольку для изопараметрического прямоугольного кольцевого
элемента Nj-1Nc = 2/5, прямоугольный элемент имеет несколько
большую степень свободы, чем треугольный. При осесимметричном
анализе использование поперечного прямоугольного элемента не
дает преимуществ. Малость отношения N, /Nc и для треугольного,
и для прямоугольного элементов не позволяет рассматривать слож-
ные деформации и одновременно удовлетворять всюду требованию
несжимаемости. Чтобы смягчить ограничение, накладываемое усло-
вием несжимаемости, необходимо вычислять скорость деформации
в одной или нескольких точках элемента.
Г. Трехмерные элементы. Тетраэдральный элемент (рис. 5.5, а)
имеет характеристики, похожие на характеристики треугольного
элемента при плоской де
нТС
рмации; компоненты скорости выражают
линейной функцией перемещения, скорость деформации постоянна
по всему элементу.
Восьмиугольный эле-
мент в виде брикета
(рис. 5.5, б) удобен,
например, при ну-
мерации элементов,
поскольку элементы
можно расположить
в определенном по-
рядке. Функцию пе-
ремещения такого эле-
Рис. 5.5. Элементы с тремя измерениями:
а — тетрагональный элемент; б - восьмиугольный блок
101
мента можно определить, разделив элемент на шесть тетраэдральных
подэлементов или используя изопараметрическую трансформацию.
5.2. Использование МКЭ для анализа
упругопластической плоской деформации
Рассмотрим треугольный элемент ijm (рис. 5.2), испытывающий
плоскую деформацию. Текущие напряжения (стх, о , а, т ); узловые
точки движутся с неизвестными скоростями (и., и.), (и , и), (ит, t>w).
Материал предполагается изотропным и подчиняется закону Гука
в упругом состоянии, условию текучести Мизеса и соотношению
Прандтля-Рейсса в пластическом состоянии. Его характеристики:
модуль Юнга Е, коэффициент Пуассона v, напряжение текучести
Напряжение текучести
и коэффициент упрочнения
Н' = dasjde
являются функциями интенсивности пластической деформации (эк-
вивалентной пластической деформацией) е = I de.
А. Скорость напряжения. Для элемента в упругом состоянии соот-
ношение между скоростью напряжения и скоростью деформации оп-
ределяется законом Гука:
(5.21)
ё, =0.
Исключая оги переворачивая вышенаписанную матрицу, полу-
чаем:
102
При пластическом деформировании полная скорость деформа-
ции {ё} включает скорость упругой деформации {ё* } и скорость пла-
стической деформации |ёр г, т.е.
(Е) = |£ J + (е J-
Упругая часть подчиняется закону Гука (5.21), а пластическая -
определяется следующими уравнениями:
£р =
Г xy J'xy 9
где ст - интенсивность напряжений (эквивалентное напряжение),
равная напряжению текучести ст в пластическом состоянии; б - ин-
тенсивность скоростей пластической деформации (эквивалентная
скорость пластической деформации).
2
(5.26)
21
£Р -£f
Р Y
ё? -£».
Согласно соотношению Прандтля-Рейсса для плоской деформа-
ции имеем следующее выражение:
— 1/2
-1/2
О'
О
О
О
ёг =0.
о
о
о
2
О
z
103
Выражения для компонентов скорости напряжений через компо-
ненты скорости деформации требуют некоторых дополнительных
уточнений, которые будут рассмотрены в разделе 5.3. Окончательное
соотношение имеет вид:
{ct}=[d₽]{4 (5.29)
Матрица
определяется уравнением
(2G)
(5.30)
где G - модуль сдвига; G = Е!2 (1 + v); , <з'у - компоненты девиатора
напряжений;
(5.31)
В общем случае соотношение между скоростью напряжения
и скоростью деформации (5.23) и (5.28) можно представить следую-
щим матричным уравнением:
{о} = [£>]{е} , (5.32)
где матрица [£>] изменяется в зависимости от образующего ее уравне-
ния.
Используя уравнения (5.32) и (5.7), скорость напряжения можно
выразить через скорости узлов:
{о} = [£>][/?]{»}- (5.33)
Б. Силы в узлах. Любое напряжение, действующее в сечении эле-
мента, должно быть уравновешено силами, приложенными вдоль бо-
ковых поверхностей. Обращаясь к рис. 5.6, находим, что составляющие
распределенной нагрузки (qml )х и (qml) на боковой поверхности т.:
(5.34)
где >5 - толщина и / - расстояние между точками т и i. Поскольку
предполагается, что силы передаются только через узловые точки, то
104
распределенные нагрузки на боковой поверхности элемента замене-
ны силами в узлах, каждая из которых равна половине равнодейст-
вующей (рис. 5.7). Составляющие силы в узле (точка i) определены
уравнениями
« = + {Цц)х^д}” (5/2){ох(у,- ~Ут)+ хху{хт — •*/)} осл
Pyi = (xm - Х})+ Хху(у; - ут )}
Силы в узлах для всего элемен-
та записывают следующим обра-
зом:
Рис. 5.6. Распределенная нагрузка
на боковой поверхности элемента
где Л =у. =у- г , ...;с~х = х -
I jnt J 7П I tfy fn
- x. ... Матрица в этом уравнении
есть транспонированная матрица
И в уравнении (5.8). Таким образом,
{р)=Л5[В]’{о),
(5-37)
где [ ]т означает транспонированную матрицу.
Рис. 5.7. Определение сил в узлах для распеделенной нагрузки (S - толщина)
105
Пренебрегая изменением формы элемента, т.е. считая А и матрицу
[В] постоянными, скорость изменения силы в узле {F} можно аппрок-
симировать уравнением
{р}= Л5[5]т{а}. (5.38)
Решая совместно (5.38) и (5.33), получаем выражение для {F} через
скорости узлов:
{р}=ладт[я][я]М (5.39)
или
И=КМ (5.40)
где [Ке]- симметричная матрица, которую называют матрицей жест-
кости элемента.
5.2.1. Набор элементов
В каждой узловой точке силы в узлах, окружающих элемент,
и внешняя нагрузка, если таковая имеется, уравновешены (рис. 5.8).
Для узла I, окруженного п элементами (/ = 1, 2,..., п) и нагруженного
внешней силой (F„.,Fy/), уравнения равновесия, выраженные в ско-
ростях изменения сил, имеют вид:
Рис. 5.8. Уравновешивание сил
и внешней нагрузки,
приложенной в узле
п
(5.41)
Поскольку скорости изменения
сил в узлах есть линейные функции
скоростей узлов (5.40), записанные
выше уравнения также являются ли-
нейными относительно скоростей
узлов окружающих элементов. На-
пример, уравнение в х-направлении
узла 7 (рис. 5.1) имеет вид:
— a|w2 +^2^2 + fl7w7 -bflgV? + +
+ a10u8 + + af2 °I1 + aBwJt2 + °14u12 - *7 -
(5.42)
106
Число N всех уравнений вдвое больше узловых точек, которые
представлены следующим уравнением, где необходимо изменить ин-
дексы
или
«1
ы2
«2/1- Fy
• ->F2, - -> fi,-. j* -+ Ao ...,Fn.
«I Fi
«11 а12 «13 a\N «2 f2
Л21 Ж a22 * «23 ж a * • a2N • ♦ b « a ► — < ►
_аЛЧ aN2 aN3 ••• aNN. • a to a a a a a
/X
[*]{«} =
(5.43)
(5.44)
где [X] - матрица жесткости.
Поскольку число линейных уравнений в системе (5.43) равно чис-
лу неизвестных параметров, система (5.44) может быть разрешена
при соответствующих граничных условиях.
5.2.2. Граничные условия
Если задана распределенная внешняя нагрузка: переднее и заднее
натяжения, силы трения и гидростатическое давление, то их исполь-
зуют для определения сосредоточенных нагрузок в узловых точках.
Увеличивающиеся скорости изменения нагрузки подставляют в пра-
вую часть уравнения (5.43). В обычных задачах формоизменения ме-
талла внутренними силами, такими как сила инерции, пренебрегают
и учитывают только действие внешней нагрузки на граничные по-
верхности.
Если заданы скорости узлов, то неизвестными считают соответст-
вующие скорости изменения нагрузок в узлах. Существуют разные
способы использования скоростных граничных условий. Поскольку
изменение в выражении для матрицы жесткости сравнительно мало,
то чаще отдают предпочтение следующему. Скорость и выражают
через С. В этом случае у-й ряд системы уравнений (5.43) приобретает
вид:
107
«11 «12 ... % — alN «1 Fl
• в • ал • • aj2 • V В • * • ... Р -аЛ в 4 • в * • ... aJN В 9 4 UJ ► =« В ♦ • р«^с/>
Л 4 4 • * « « • « в • • 4 « » • V в • в • 4
«М aN2 - «^ — aNN .UN. Fn
(5.45)
где р - очень большое число (108...1012).
Из j-го уравнения (5.45) следует, что и. = С. Решив уравнение,
можно определить скорость изменения силы F в узле:
Ру=Уа}кик. (5.46)
fc=l
Если в случае, представленном на рис. 5.9, предполагается нали-
чие коэффициента трения ц, то к уравнению (5.43) следует добавить
следующие два уравнения:
иА-А
= (- ц sin 6 + cos 6)FX+(- pcos 0 — sin 6)Fy
и
и = и sin 0 - о cos 6 = 0.
и
(5.47)
Рис. 5.9. Кинематика и силы при течении мегалла вдоль наклонной
плоскости с заданным коэффициентом трения
108
5.2.3. Решение матричного уравнения для приращений
Приведенные выше линейные уравнения решаются с использова-
нием ЭВМ одним из методов типа процедуры подметания Гаусса или
процедуры Гаусса-Зейделя. Члены матрицы [А] отличны от нуля
обычно только вблизи главной диагонали, и, следовательно, при ре-
шении большого матричного уравнения на ЭВМ сохраняются только
члены, собранные вблизи диагонали.
После вычисления скоростей узлов согласно уравнениям (5.7)
и (5.33) можно вычислить скорости деформаций и скорости измене-
ния напряжений. Компоненты скорости пластической деформации
находят из соотношения
Скорость эквивалентной пластической деформации задается
уравнением (5.27). Приращения определяют как произведения скоро-
сти на интервал времени А/ и добавляют к предыдущим значениям,
находя, тем самым, новые значения напряжений, координаты и ин-
тенсивности (эквивалентные) пластической деформации. Эти новые
значения затем используют как исходные при вычислениях на сле-
дующем шаге.
5.2.4. Упругопластический переход
В процессе нагружения на разных стадиях может начинаться пла-
стическое течение материала. Если приращение деформации или
временной интервал оставить постоянными, то после малого прира-
щения деформации напряжение в упругом элементе может превысить
напряжение текучести. Чтобы избежать этого превышения, времен-
ной интервал следует так модифицировать, чтобы интенсивность
напряжений в элементе совпадала с напряжением текучести после
малого приращения деформации. Когда число элементов велико,
процедура требует очень длительного времени для проведения вы-
числений, чтобы текучесть в элементах наступала последовательно.
Для сокращения времени вычисления упругие элементы, интенсив-
ность напряжений в которых после очередного приращения близка
к напряжению текучести (99,5 % или более), будем рассматривать как
уже достигшие текучести на этом шаге.
При разгрузке или прекращении в некоторых элементах пластиче-
ских деформаций элементы вновь деформируются упруго. Чтобы
обнаружить существование такого перехода, для каждого элемента
109
вычисляют скорость диссипации энергии пластического деформиро-
вания:
w.
z^z
Если диссипация энергии оказывается отрицательной (5.49), то
считают, что элемент вернулся в упругое состояние, и [Z)j заменяют
на
5.3. Теория МКЭ для упругопластической
деформации
5.3.1. Получение матрицы [ор]
А. Обобщенная функция текучести. Для получения матрицы [£>р]
необходимо провести инверсию соотношения между напряжением
и деформацией.
При пластической деформации упругопластического материала
скорость деформации {ё} состоит из упругой и пластической частей:
ёе
В общем случае упругая часть подчиняется закону Гука:
где точное выражение [/>* ] для изотропного материала дается для
трехмерной деформации в декартовых координатах
щим выражением:
f(l-v)/(l-2v)
z следую-
v/(l-2v)
О
О
О
v/(l-2v)
(l-v)/(l-2v)
v/(l-2v)
О
О
О
v/(l-2v)
v/(l-2v)
(l-v)/(l-2v)
О
О
О
о
о
о
1/2
О
О
о
о
о
о
1/2
О
О
О
О
О
о
1/2
Считая, что пластический потенциал материала равен /, пластиче-
скую часть можно записать так:
^}=л{а//ао}/, (5.53)
ПО
где h (>0) - константа. В дальнейшем предполагаем, что пластиче-
ский потенциал можно отождествить с функцией текучести, что явля-
ется общим случаем:
/=о, (5.54)
в что он монотонно увеличивается (упрочняется) с увеличением ин-
тенсивности (эквивалентной) деформации ё. Прежде чем выводить
матрицу ], запишем некоторые соотношения между Л {Э/-/да} и /.
Введем модуль упрочнения И*, связанный с f соотношением
f = b = H't. (5.55)
Подставляя (5.55) в (5.53), получим
{ё'}=й{а//да}Н'1. (5.564
Мощность диссипации определяется выражением
Wp ={о}т {ё/’}=стё = Л{а}т{5//5а}Ягё, (5.57)
и, следовательно,
\lh = tya){a}T{df/da}H'. (5.58)
Запишем полный дифференциал от/:
df = (df/дах}дах + (df/day^ay + ...+(5//^2)^г + ... =
= (df / da)T да, (5.59)
или, выражая через скорости изменения напряжений, получим
f -{df /Эо}{ст}.
(5.60)
Теперь можно выразить (о) через {ё}. Объединяя уравнения (5.51)
и (5.53), получим
{о}ф* ]({ё}- {ё' })= [/Г ] ({ё}- h{df/da}f\ (5.61)
Подставляем это уравнение в (5.60):
/ = {df/da}r [ое ] {ё* }- h{df /да}т [ое ] [df / ост}/. (5.62)
111
Таким образом,
А/ {У/а°)т к ] И
-+{§г/аа}т к ] {а//Эа}
h
Уравнения (5.53) и (5.61) совместно с (5.63) дают
Ы = ИМ’
гр, 1 rDe 1 к к//Мд//Мт к ]
(1/ст){а}т {д//ди}Н’ + {df/do}' к ] {df/da}'
Б. Условие пластичности Мизеса. Рассмотрим изотропный упроч-
няемый материал, подчиняющийся условию пластичности Мизеса:
/2=п2={(о>,-ог)2 +(ct.-ctx)2(ctx-ctj,)2 +6(tJz+t2x+t2,)}/2. (5.65)
Из этого уравнения получаем
{а//3о}т = з{стхо^ст'г2тгг2тгх2тх?}/2ст,
где а'х, ст' , ст. - компоненты девиатора напряжений.
Используя уравнение (5.52), получаем
С</ / j y^z yz ^zx хуг
где G- модуль упругости второго рода, G = E/ 2(1+v).
{df /аст}т к ]{3/7 аст}=3G,
к ]{df / dc}{df / ди}т к ]=
f
9G2
(5.66)
(5.67)
(5.68)
(5.69)
(5.70)
Z
yz
zx
112
Таким образом,
М= [D' ]________________
4о2(ЗС + Я')/9
В. Матрица [ZX] для некоторых типов деформаций. Для плоской
деформации ё. = y2X = yyz = О и т * = т = 0. Матрица [D ] имеет 3x3
компонент, соответствующих стх, Ьу, тх>, ёх, ъу, уху, и получается
путем исключения ненужных компонент из матрицы 6x6, а именно
из уравнений (5.52) и (5.71). стг можно также вычислить, пользуясь
выражениями для ёх, ё.,, уж?.
Для плоского напряженного состояния ст2 = iyz = izx =0
и у „г = у2Х = 0. При получении матрицы [яе ] необходимо выделить д2 .
В результате этой процедуры компоненты матрицы будут отличны
от тех, которые получаются из уравнения (5.52). При плоском на-
пряженном состоянии матрица [De I для условия пластичности Мизе-
са имеет вид
v/(l + v)
1/(1 ч-v)
0
0
0
(l-v)/2(l + v)
1/Cl + v)
D'l=----- v/(l + v)
i-v 0
® х^ху
У ХУ
(5.72)
Для осесимметричной деформации координаты х, у, z заменяют
на г, в, г. Если отсутствует скручивание, что типично для большинст-
113
ва задач деформирования металлов, то = увг = 0 и = т0г = 0.
Матрица [/)] содержит 4x4 компонент, соответствующих сгг,дв,
®г»^гт>^г» ^6’ ®г» Угг-
5.3.2. Определение силы в узле методом виртуальной
работы
Рассмотрим произвольный набор скоростей узлов ) элемента
и соответствующее распределение скоростей деформации внутренне-
го {ё* г. Принцип виртуальной работы выражается уравнением
|{ё‘){а}<7Г = {к‘}{р},
(5.73)
где {о} и {р} - действительное напряжение и реальная внешняя на-
грузка соответственно.
Это уравнение означает, что скорость диссипации внутренней
энергии на виртуальных деформациях равна мощности внешних сил
на виртуальных перемещениях. Подставляя соотношение
ИИ’И’
(5.74)
в уравнение (5.73), получаем
k}Tf !«’)Т{р}; (5-75)
{p} = J[2?]T{a}JK. (5.76)
Если изменением формы элемента, т.е. изменением в матрице [В]
и изменением объема пренебрегают, то скорость изменения силы
в узле [р] задается соотношением
(5.77)
После подстановки в уравнение (5.77) соотношения
{б}= №}
(5.78)
оно приобретает вид:
114
W=J [вГ к»)' - {Jw И№и} М=
-{|[*ки}{в}=к,]{«}.
(5.79)
Если матрица [В] неоднородна внутри элемента, то интегрирова-
ние не всегда удается выполнить аналитически. Практически выпол-
няют численное интегрирование путем выбора нескольких Гауссовых
точек в элементе с соответствующим выбором их «веса». На рис. 5.12
дан пример выбора Гауссовых точек и их «веса» для треугольника
и тетраэдра. Матрица [Л?е ] для элемента объемом V аппроксимиру-
ется «-Гауссовыми точками следующим образом:
= (5.80)
Простейший способ - оценивать матрицу [J3] только в центре. На-
пример, осесимметричного элемента можно аппроксимировать
так:
M = 2K[8]T[Z)][s]r4, (5.81)
где г и В - соответственно радиус и значение [В] в центральной
точке; Ас - площадь поперечного сечения.
Рис. 5.12. Точки Гаусса и их «вес»
Известно, что аппроксимация дает удовлетворительный резуль-
тат, если число элементов достаточно велико.
5.3.3. Получение матрицы жесткости с использованием
вариационного метода
В предыдущих разделах основополагающие уравнения типа (5.44)
были получены из условия равновесия в узловых точках. Те же урав-
115
нения можно получить с помощью вариационного метода, не обра-
щаясь к понятию силы в узле.
Любое поле скоростей, которое удовлетворяет уравнению совмес-
тимости (соотношению между скоростью и скоростью деформации)
и скоростному граничному условию, называется кинематически воз-
можным. Рассмотрим тело V, поверхность которого состоит из час-
тей Su и S: скорости определены на Su, поверхностная сила сцепления
{<?} определена на S. Если выбирается кинематически возможное
поле скоростей {и}, то соответствующая скорость изменения напря-
жений
{о}=[/)]{ё} (5.82)
не всегда удовлетворяет граничному условию для напряжений и ус-
ловию равновесия, она может не оказаться статически возможной.
Первый вариационный принцип гласит: среди кинематически
возможных решений точное решение соответствует абсолютному
минимуму функционала
J3=O,5j{a}T |{tf}T {v}dS.
? St
(5.83)
Этот принцип предполагает, что минимизация функционала по-
зволяет определить точное решение. Модель конечного элемента
можно использовать для нахождения минимума функционала, хотя
это решение может оказаться не совсем точным из-за дискретности
представления.
Поскольку {ё} и {<у} - линейные функции скоростей узлов и, то
{ё} = [В] {«}, (5.84)
(5.85)
и
{о}т = = {м}т[В]т[Р]. (5.86)
Находя равнодействующие распределенной поверхностной на-
грузки {q} в узловых точках (узловая сила {F}), можно выразить
функционал уравнения (5.83) через скорости узлов:
узел на
]J= £о,5{а}т
элемент
{к}. (5.87)
116
Минимизацию можно выполнить, дифференцируя функционал по
скоростям узлов каждого элемента. Легко понять, что вычисления
упрощаются из-за симметричного характера матрицы [Ке ]. Условия
стационарности (минимума) записываются так:
И4-°- (5.88)
Это уравнение тождественно уравнению (5.44).
5.3.4. Теория конечной деформации
В предыдущих разделах силу в узле после бесконечно малой де-
формации определяли по теории, в которой предполагалось, что
форма или матрица [В] каждого элемента не изменяются при беско-
нечно малой деформации. Это называется инфинитезимальной тео-
рией деформации. Для получения более точного результата необхо-
димо рассмотреть влияние изменения формы и вращения элемента по
времени каждого приращения малой деформации. Такая теория на-
зывается теорией конечной деформации. В дальнейшем теорию ко-
нечной деформации рассмотрим на примере треугольного элемента
при плоской деформации (площадь А и толщина S).
Поскольку силы в узлах заданы уравнением (5.37) в виде
{ст}, то полное выражение для скорости изменения силы
в узле
{/>} = Л5[В]ТН (5-89)
Скорость изменения напряжения {ст} в этом выражении зависит не
только от деформации (5.29), но также и от вращения элемента. При
вращении элемента, даже если напряжение, связанное с материально
фиксированными координатами, остается постоянным, напряжение
Эйлера или Коши {ст}, связанное с пространственно фиксированны-
ми координатами, изменяется. Очевидно, скорость изменения напря-
жения Йомена р7 ’, используемая для характеристического уравне-
ния, должна определяться материально фиксированными координа-
тами. Таким образом,
р}=№}- (5.90)
В случае двумерной деформации {ст} и 4г7 } связаны соотношени-
ем, содержащим угловую скорость вращения элемента против часо-
вой стрелки со:
117
(5.91)
Таким образом, уравнение (5.89) приобретает вид
{/>}= т{о) + Лв[в]т{о}+ ' }+<оху
(5.92)
При пластическом деформировании обычных металлов влиянием
изменения объема можно пренебречь:
{/>} = A s[b]t {о} + аху Л5[В]
ЛЯ[В]Т[П]№). (5.93)
Сравнивая это уравнение с уравнением (5.39), можно увидеть, что
инфинитезимальная теория деформаций учитывает только последний
член в этом выражении.
Поскольку [в] и со.^ являются функциями узловых скоростей, со-
гласно уравнениям (5.8) и (5.13), то уравнение (5.93) имеет следую-
щий вид:
{/>} = К«} + [К2 &} + Ы4 (5.94)
где [Kj] - матрица, связанная с изменением формы элемента; [АГ2] -
матрица, связанная с вращением; [АГ0] - матрица, связанная с измене-
нием деформации и идентичная матрице [X] из уравнения (5.44).
Явный вид [ЛГ^и |Х2]:
О 0 - т,„
АЛ
О 0 -<уу
ст т
^х ху
^ху
О ~Тху
о -оу
О
~*ху О
(5.95)
118
и
где
3
8Л
~e\ci
c6^i cfS-j
c^j — схст
^2^j ~ ^2^т
^3^J ~~ ^З^т
C4 bj C4Cm
C$bj ~ c5cm
C(>bj ~ c6cm
c4^m
(5.96)
(5.97)
м=
Это выражение, учитывающее изменение матрицы [3], позволяет
выразить силы в узлах после деформации через координаты после
деформации.
5.4. Использование МКЭ для анализа
жесткопластических деформаций
5.4.1. Общие сведения
Для упрощения рассуждений при анализе деформирования металл
считают жесткопластическим, поскольку упругие деформации, как
правило, намного меньше пластических. При таком предположении
напряжения в жесткой зоне определить нельзя. Более того, для не-
сжимаемых материалов компоненты гидростатического давления
также не связаны с компонентами скорости деформации. По этой
причине анализ деформирования жесткопластических материалов
МКЭ отличается от анализа упругопластических.
5.4.2. Вариационный принцип Маркова
Рассмотрим тело из жесткопластического материала, удовлетво-
ряющее условию пластичности Мизеса (уравнение 5.59):
119
Соответствующее правило течения пластической деформации оп-
ределяется уравнениями Леви-Мизеса, в котором скорость деформа-
ции {ё} связана с компонентами девиатора напряжений {ст'}, т.е.
рмации.
(5.98)
где е - интенсивность скорости де
Уравнение (5.98) удовлетворяет условию несжимаемости
ёг = ё„ + ё„ +ё, =0.
г -Л Jr
(5.100)
Напряжения выражаются через скорость деформации с помощью
гидростатического давления:
= о7+а™
СТ- = ст'- Ч-Ст —
z z т
(5.101)
Поскольку гидростатическое давление не связано со скоростью
деформации, то нормальное напряжение нельзя определить одно-
значно по скорости деформации, кроме случая, когда гидростатиче-
ское давление или одно из нормальных напряжений являются задан-
ными, как это бывает при плоском напряженном состоянии. Эта воз-
можность дается условием несжимаемости.
Предположим, что объем тела V, а площадь поверхности S состо-
ит из Su и S частей: скорости определены на и внешние силы {q} -
120
на S . Любое поле скоростей, удовлетворяющее условию несжимае-
мости в теле и граничным условиям в скоростях на границе Su, явля-
ется кинематически возможным. Вариационный принцип Маркова
гласит: «Среди всех кинематически допустимых полей скоростей,
точное решение обеспечивает абсолютный минимум функционала
Ф = [ aedF - J {?}т {u}dS ». (5.102)
И S,
Этот принцип является альтернативным выражением теоремы
о верхней оценке. Первый член функционала (5.102) есть скорость
рассеяния (диссипации) энергии пластической деформации, второй
член - работа, совершаемая внешними силами сопротивления. Если
отдельно определено касательное напряжение, вызываемое трением,
как в большинстве задач деформирования металла, то функционал
состоит из диссипации энергии пластической деформации и потерь
на трение.
Из этого принципа следует, что точное поле скоростей можно оп-
ределить путем минимизации функционала Ф. Использование конеч-
ных элементов целесообразно для нахождения достаточно прибли-
женного поля скоростей. Пользуясь уравнением (5.7), интенсивность
(эквивалентную) скорости деформации ё можно выразить через ли-
нейную функцию скоростей узлов. Первый и второй члены уравне-
ния (5.102) определяются соответственно так:
Ф1 = e(ui,u2, ...,u„)dy (5.103)
элемент
И
Ф2=- £{F}t{4 (5.104)
узел на 5,
Подставляя эти соотношения в (5.102), получаем
Ф = ё(и!,н2,..., u„)dV-£{^}т{ы}. (5.105)
Условие несжимаемости каждого элемента также можно выразить
через скорости узлов:
V = EyV = У\ик = 0, (5.106)
где Ёу - скорость средней деформации изменения объемау-го элемента.
121
Для треугольного элемента при плоской деформации (рис. 5.2) яв-
ный вид уравнений (5.105) и (5.106) можно получить, пользуясь соот-
ветственно уравнениями (5.8) и (5.10).
+ 3(<од + + CyUj + b^j + + bmvm)2/1}/9 - {F)T{«}. (5.107)
(5.108)
Задача оптимизации в методе верхней оценки с помощью МКЭ зву-
чит так: найти набор скоростей узлов, минимизирующий нелиней-
ный вещественный функционал уравнения (5.105) при заданных ли-
нейных ограничениях уравнения (5.106) и граничных условиях
в скоростях.
5.4.3. Методы оптимизации
А. Метод линейного программирования. Метод линейного про-
граммирования создан для определения максимального или мини-
мального значения линейной вещественной функции при линейных
ограничениях.
Вещественная функция G = для которой определяют мак-
симум или минимум при ограничивающих условиях:
0 = 1>2,...,ли2>
<0 U = 1> 2> —."«з),
(5.109)
где х. (i = 1, 2, ..., т) - положительные переменные. В практических
целях для решения задач линейного программирования используют
такие методы, как метод симплекса и дуплекса.
Поскольку вещественная функция уравнения (5.104) нелинейная,
необходимо линеаризовать ее с помощью затравочного поля скоро-
стей, чтобы в дальнейшем оперировать методом линейного про-
граммирования. Предполагая, что набор скоростей узлов, аппрокси-
мирующий точное решение, есть ы01, и02, ..., можно получить опти-
122
мальное решение, добавляя к ним небольшие скоростные поправки
dult du2 и т.д. Линеаризованные вещественная функция и ограничения
имеют вид:
Ф = Ф0+(дФ/ди1)0^1 +(дФ/ди2)0<1и2 +..., (5.110)
и
где (дФ/диу )0, (дФ/ди2)ц и Фо - значения, вычисленные с помощью
затравочного поля скоростей.
К затравочным значениям прибавляют поправки, полученные ме-
тодом линейного программирования. Новые значения скоростей,
в свою очередь, используются в качестве затравочных на следующем
шаге итерации. Эта процедура продолжается до тех пор, пока по-
правки не обратятся в нуль.
Б. Метод дополнительной функции (пенальти). Большое положи-
тельное число р умножаем на квадрат скорости изменения объема
и прибавляем к выражению функционала (5.105). В результате полу-
чаем новый функционал
у=Ф + р22(ёгг)2. (5.111)
i
Если в каком-нибудь из элементов отлично от нуля, то приме-
няют метод дополнительной функции (пенальти). Когда новый
функционал достигает минимума, скорость деформации изменения
объема будет близка к нулю.
Поскольку новый функционал у свободен от разных ограниче-
ний, его минимальное значение проще найти, чем решить первона-
чальную задачу. Простым, но отнимающим много времени методом
является метод последовательного приближения, т.е. изменение ско-
ростей узлов в направлении, соответствующем меньшему значению
функционала. Более изящным, возможно, является следующий метод.
Дифференцируем функционал по скоростям узлов и результат при-
равниваем к нулю:
ду/дщ =/i(«t,H2, ...,ил) = 0;
ду!ди2 =/2(щ,и2, ...,м„) = 0, (5.112)
Нелинейные уравнения можно решать методом итераций, кото-
рый называют методом Ньютона-Рафсона. С помощью затравочных
скоростей и поправок к ним линеаризуют уравнения для f^f2>
123
A = U)o + VW2 + ->
/2 = (/2 )0+(ЭЛ fax V“i+ (36 /aw2 )0Ж2 + • •,
(5.113)
где /^2)0’— - величины, вычислен-
ные с помощью затравочного поля скоростей. Решая линеаризован-
ные уравнения, определяем поправки для скорости. Используя новое
поле скоростей в качестве затравочного, продолжаем итерацию, пока
результат не сойдется.
5.4.4. Вычисление напряжений в очаге деформации
А. Метод интегрирования. Изучение метода визиопластичности
показало, что напряжение в очаге деформации можно вычислить
с помощью интегрирования дифференциального уравнения. В случае
плоской деформации следующие уравнения получаются из условия
равновесия:
у
= - f &Ху !дх + дъ'у /ду)х^хо dy + ото ,
Уо (5.114)
«т = - f /дх + да'х /ду)у=уо dy + Ото ,
*о
где ст_ - известное значение в точке х = х0, у =у0.
При применении этого метода в МКЭ используют компоненты
девиатора напряжений, связанные с оптимальным полем скоростей.
Этот метод, как правило, дает большой разброс значений напряже-
ний, поскольку произвольное деление на конечные элементы не все-
гда позволяет точно вычислить градиенты компонент девиатора на-
пряжений бах/йх а поэтому интегрирование по разным путям дает
разные значения для напряжений.
Б. Метод множителей Лагранжа. Метод множителей Лагранжа -
это классический подход к задаче оптимизации нелинейной вещест-
венной функции при наличии ограничивающих условий. Принцип
этого метода заключается в следующем: «Оптимизация вещественной
функции б(х1,х2,...,хи)-^ максимум или минимум при наличии ог-
раничивающих условий
g, (х,,х2, ...,хл)=0 (i=l,2, ...,m) (5.115)
124
эквивалентна определению условия стационарности функционала
О — G(xj , Х2 * хл )+ X j (^1 > ^2 > —» )» О- 1
i
где X (i = 1,2,..., т) - множители Лагранжа».
Условие стационарности определяется следующими уравнениями:
д£1/дщ=О; дС1/ди2-О; ... дС1/ди„ =0; fS 1171
SQ/SX^O; dQ/dX2 =0; ... Xl/&km =0. Р‘ }
п + т неизвестных параметров хр х2, .... х и Х(, ..., X определя-
ются в результате решения п + т уравнений. При этом т уравнений
являются частными производными по Хь Х2, ..., X и совпадают с ог-
раничивающими условиями. Таким образом, члены X^j, X^g, ...
в уравнении (5.116) обращаются в нуль, когда функционал уравнения
(5.116) достигнет минимума. Отсюда заключаем, что минимальное
значение О при данных ограничениях равно минимальному значе-
нию G.
В рассматриваемой проблеме вещественной функцией является
уравнение (5.105), а ограничением - условие несжимаемости (5.106).
Функционал О в нашей задаче имеет вид:
(5.118)
Условия стационарности
/1(и1>и2» ”-»Х1>
f =F
J п л n
(5.119)
и Vi(«i,u2> • •)=0,v2 = 0,om =0.
Уравнение (5.119) решается методом Ньютона-Рафсона с помо-
щью затравочного поля скоростей и их поправок. Характерная черта
этого метода - равенство множителей Лагранжа компонентам гидро-
статического давления при достижении функционалом минимума.
Так вычисляют напряжения. Это можно показать следующим об-
разом.
Предположим, что этим методом определяют точное поле скоро-
стей. Функционал Ф, выраженный через напряжение {о} и скорости
деформации {ё}, приобретает вид:
125
ф=X f {°F - X w=X J ^dv+
i J *
+Xf°'"^r
I
“X W = X f ° ® dV -
j i
" Х^ + X J °’
j I
(5.120)
где {o'} обозначает компоненту девиатора напряжений, а - гидро-
статическое давление (среднее напряжение).
Если гидростатическое давление предполагается постоянным
в каждом элементе, то уравнение (5.120) записывают так:
Ф = S j И’ {“} + X °™!«i
i J •
(5.121)
Сравнивая это уравнение с уравнением (5.118), видим, что множи-
тель Лагранжа представляет гидростатическое давление в том случае,
если полученный результат точен. Таким образом,
^ч' °пи •
(5.122)
Другое доказательство возможно из рассмотрения равновесия сил
в узлах. Силу в узле с помощью уравнения (5.76) можно записать
{Р} = J [Бр {a}dV = J [Б]7 {о ’}dV• + J [б]
'т
0
О
о
МИ.
(5.123)
Поскольку компоненты девиатора напряжений выражаются через
скорости узлов (5.102), то силы в узлах элемента и неизвестное гидроста-
тическое давление можно записать через скорости узлов. Уравнения рав-
новесия сил в узловых точках имеют вид /](«1,Ы2» •••> СТЛ11» СТт2> •••) = Л’
Уз ~ Уг»• • •»Jn ~ Fn •
Эти уравнения аналогичны (5.119) при fx — Гх,/г = F2, — Fn ,
если Xj, ... заменить на стт1, от2» •••
126
В. Метод сжимаемости. Рассмотрим условие пластичности в форме
о* =^0^|оу-аг)г +(ог-стх)2 +(стЛ-о>,)2 + б(т^г + + т^)}+ао£,. (5.124)
Условие пластичности (5.124) учитывает влияние гидростатиче-
ского давления. При а, равном нулю, уравнение (5.124) совпадает
с условием пластичности Мизеса. Чтобы аппроксимировать условие
пластичности Мизеса, значение а должно быть очень мало.
Учитывая соответствующее правило течения (5.124), напряжения
можно выразить через деформации:
где е* — интенсивность (эквивалентная) скоростей деформации:
(5.126)
Из-за чувствительности условия пластичности к гидростатиче-
скому давлению материал при пластической деформации претерпе-
вает изменение объема. Благодаря этому можно вычислить непосред-
ственно напряжения, зная скорости деформации.
В табл. 5.3 показано относительное отклонение модифициро-
ванного условия пластичности от условия пластичности Мизеса при
некоторых значениях а. Отклонение напряжения текучести при про-
стом сжатии определяется членом (о* -ст)/а, а объемные деформа-
ции оцениваются при простом сжатии до степени деформации 10 %.
При а меньше 0,01 эти отклонения очень малы.
Таблица 5.3
а * «в» ст — ст а Изменение объема, %, после осадки на 10 %
0,04 0,2 -0,1
0,01 0,06 -0,04
0,0025 0,01 -0,009
127
Вариационный принцип записывают в таком же виде, как и урав-
нение (5.102), за исключением того, что интенсивность напряжений
и интенсивность скоростей деформации включают соответственно
гидростатическое давление и скорость деформации объема:
Ф* = f о *ё W - f {?}т {u}d$. (5.127)
г s,
Условие несжимаемости не является необходимым условием для
кинематически возможного поля скоростей, и, следовательно, миними-
зацию функционала (5.127) можно провести, не используя условия по-
стоянства объема. Условие стационарности получается просто путем
дифференцирования вещественной функции по скоростям узлов:
ЭФ7Эи2=0, 5Ф7ам„=0. (5.128)
Окончательный результат получаем методом Ньютона-Рафсона.
Объем каждого элемента сохраняется почти постоянным, поскольку
изменение объема ёи ведет к большому увеличению интенсивности
скоростей деформации и диссипации энергии. Член, соответствую-
щий изменению объема, играет ту же роль, что и дополнительный
член в методе дополнительной функции. Преимущество этого метода
в том, что анализ с использованием МКЭ в значительной степени
упрощается в результате возможности непосредственного вычисле-
ния напряжений по скоростям деформации.
5.4.5. Некоторые аспекты использования МКЭ
для жесткопластического материала
А. Затравочное поле скоростей. Затравочное поле скоростей необ-
ходимо для процедуры оптимизации при использовании МКЭ для
жесткопластического материала. Из опыта известно, что выбор за-
травочного поля скоростей существенно влияет на продол-
жительность счета. Неподходящее затравочное поле скоростей во-
обще может не дать сходящегося результата.
Затравочное поле скоростей можно получить простым способом,
например в случае простой геометрии заготовки методом верхней
оценки. Однако в общем случае выбор правильного затравочного
поля скоростей для практических проблем - задача непростая, и не-
обходима процедура, позволяющая находить подходящее затравоч-
ное поле скоростей.
Рассмотрим функцию £?, похожую на функцию в уравнении
(5.118),
128
б = ^(о„ё1,Г<)2 +22±(f7u7J! +£а.Х-_
(5.129)
Знак «+» соответствует FjVj 0 и «-» FjVj < 0. Минимизировать Q
достаточно просто, поскольку частные производные по скоростям узлов
и 1 являются линейными. Это решение обычно дает достаточную
точность для оптимального поля скоростей.
Б. Критерии сходимости. Процедура нелинейной оптимизации
требует итерационных вычислений для получения оптимального ре-
шения из затравочного поля скоростей. Итерации обычно продол-
жают до тех пор, пока поправки скорости не станут очень малы, но
этот метод не дает гарантии, что полученный результат обеспечит
хорошую точность. Поскольку всегда выбирают кинематически воз-
можное поле скоростей, решение будет верным, если соответствую-
щее поле напряжения является статически возможным. Из этих сооб-
ражений предложен критерий сходимости, основанный на равнове-
сии сил в узлах.
Если установлены напряжения в элементе, то силы в узлах опре-
деляют из уравнения (5.76). Из рис. 5.14 видно, что неуравновешен-
ная часть сил в узловой точке i при двухмерной деформации
(5.130)
Можно оценить степень несба-
лансированности во всей облас-
ти AS:
где N„- полное число узловых
точек и А - среднее значение
Рис. 5.14. Силы в узлах элемента
площади граничных поверхно-
стей элемента.
Если величина AS очень мала
по сравнению с напряжением те-
кучести, 1/100 или меньше, то можно считать, что решение для поля
скоростей обеспечивает по отношению к оптимальному хорошую
точность.
129
В. Выбор элемента. Как указывалось ранее, ограничение, накла-
дываемое условием несжимаемости, уменьшает степень свободы де-
формации.
Рассмотрим случай, когда на углу заготовки при осадке выбира-
ют треугольный элемент (рис. 5.15, а). Нельзя получить значение уг-
ла, увеличивающееся при повороте боковой поверхности от 90°
до 180°, без уменьшения объема углового элемента до нуля. Четырех-
угольный элемент (рис. 5.3, б) имеет большую степень свободы, чем
треугольный, если условие несжимаемости рассматривается в точке
каждого элемента, поэтому его можно использовать на угловых
кромках заготовки (рис. 5.15, б). В практических расчетах при вы-
числении уравнений (5.99) и (5.127) используют скорость деформации
объема в центре элемента (^ = 0, ц = 0) независимо от положения
элемента, в то время как другие компоненты скорости деформации
можно варьировать внутри элемента.
Г. Конечные дефор-
мации и упрочнение.
В рассмотренной выше
формулировке МКЭ для
жесткопластического ма-
териала влияние упроч-
нения учитывалось толь-
ко для распределения
потока напряжений. Ма-
териал в процессе при-
ращения деформации
(рис. 5.16) предполагали
идеально пластическим,
и для учета упрочнения
в конце каждого шага
использовали экстрапо-
ляцию в начале шага из рассмотрения функционала или условия рав-
новесия сил в узлах. Однако в некоторых условиях деформирования
деформация оказывается неустойчивой из-за образования шейки,
потери устойчивости и появления локального сдвига. Таким обра-
зом, на полученный результат оказывают влияние изменение объема
и характеристики упрочнения. При рассмотрении таких задач необ-
ходимо учитывать конечную деформацию и коэффициент упрочне-
ния на шаге деформации.
а б
Рис. 5. /5. Угловые элементы при осадке заго-
товки (/ — инструмент): а - треугольный элемент;
б - четырехугольный элемент
130
Рис, 5,16. Схема учета упрочнения МКЭ
для жесткопластической деформации: а -
упрочнение отсутствует на шаге; б - упрочнение,
линейное на шаге
Один из возможных методов - это составление уравнения равно-
весия в конце каждого приращения де
Hit
рмации путем рассмотрения
изменений форм элементов и коэффициентов упрочнения (рис. 5.16, б).
При решении этой задачи методом минимизации функционала сле-
дует использовать напряжения, определенные для больших деформа-
ций, т.е. напряжения Кирхгоффа.
5.5. Применение МКЭ в процессах штамповки
(вопросы, возникающие при использовании МКЭ)
А. Характеристическое уравнение. При анализе процессов штам-
повки металлы часто предполагают однородными и нечувствитель-
ными к скорости, удовлетворяющими условию пластичности Мизеса
и связанному с ним правилу текучести, хотя реальные металлы всегда
обнаруживают некоторые отклонения от этих предположений. Если
отклонение становится значительным, то следует соответствующим
образом изменить характеристическое уравнение.
Чувствительность к скорости рассматривалась при анализе сверх-
пластического материала. При высоких температурах металл можно
считать вязким или вязкопластическим. В рамках жесткопластиче-
ского материала МКЭ легко учесть влияние скорости, в то время как
в рамках упругопластического материала для этого требуется итера-
ционная процедура.
Б. Граничные условия с трением. Касательные напряжения меж-
ду поверхностями образца и инструмента определяются многими
факторами: сочетание металлов, скорость скольжения, смазочный
материал и др. Точное выражение, описывающее трение, если оно
и существует, крайне сложно. Однако в практических целях можно
использовать простые зависимости, учитывающие влияние трения.
Для этого обычно используют коэффициент трения или фактор
трения. При анализе МКЭ можно применить правило скольжения,
131
предложенное Сегучи и др. В этом случае удобно рассматривать
скольжение тонких слоев с учетом трения.
В. Сингулярные (особые) точки. В большинстве технологических
задач штамповки существуют области, в которых направление тече-
ния меняется резко, как показано на рис. 5.17, а (точки А и Б). Если
эти области рассматривать как сингулярные (особые) точки, то, как
правило, МКЭ, применяемый для анализа штамповки, будет трудо-
емким. Один из способов преодоления этих трудностей - рассматри-
вать очень тонкие элементы вокруг точки. Однако это требует увели-
чения продолжительности счета. Для получения разумно хорошего
результата без увеличения числа элементов можно использовать эле-
мент, как показано на рис. 5.17, б.
Рис. 5.17. Сингулярные (особые) точки:
а - схема течения: 1 - матрица, 2 - контейнер; б - четырехугольный элемент; в -
треугольный элемент
В качестве экстремального случая двух близко расположенных уз-
ловых точек, одна - на поверхности матрицы, другая - на поверхности
контейнера, которые движутся в разных направлениях, рассматривает-
ся узел, обладающий двумя скоростями. Элемент с сингулярной
(особой) точкой можно рассматривать так же, как изопараметриче-
ский четырехугольный элемент. Результаты расчетов изменения объе-
ма и напряжений улучшаются, если в МКЭ для анализа жесткопласти-
ческой д
рмации использовать такой элемент.
Г. Стационарные процессы. Для анализа непрерывных процессов
типа выдавливания, вытяжки и прокатки с использованием непре-
рывных итерационных вычислений от начального до стационарного
состояния требуется очень большая продолжительность для прове-
дения вычислений. Решение, соответствующее стационарному со-
стоянию, можно получить также путем модифицирования линий то-
ка. Предложенные вначале линии тока модифицируют с помощью
результирующего поля скоростей. Эту модификацию продолжают до
тех пор, пока поле скоростей не сойдется. Значение напряжения теку-
чести в каждом элементе получают путем интегрирования интенсив-
ности скорости деформации вдоль линий тока элемента. Эта проце-
132
дура особенно целесообразна при использовании МКЭ для анализа
жесткопластических, а также упругопластических деформаций.
Анализ жесткопластических деформаций МКЭ основывается на
методе множителей Лагранжа. Абсолютный уровень гидростатиче-
ского давления неопределим в случае, когда зоны пластичности ок-
ружены жесткими зонами и штампом, а граничные условия для на-
пряжений неизвестны. Во избежание этого используют специальный
прием, а именно, уровень гидростатического давления определяют из
условия равенства полной силы на входе или выходе и внешней силы
(заднее или переднее натяжение). Метод, основанный на свойстве
сжимаемости материала, при рассмотрении непрерывных процессов
с такими трудностями не сталкивается.
Д. Продолжительность счета. Продолжительность счета значи-
тельно уменьшилась в результате увеличения быстродействия ЭВМ.
Однако скорости все же недостаточны для практического примене-
ния метода к процессам штамповки, поскольку итерационное вычис-
ление больших матричных уравнений приходится повторять сотни
раз, чтобы описать конечную деформацию. Уменьшение числа эле-
ментов и увеличение размера каждого шага существенно уменьшают
продолжительность расчетов, хотя это может привести к ухудшению
точности решения, особенно в случае упругопластического анализа.
Для анализа жесткопластического материала можно применять
сравнительно большой шаг, поскольку значения напряжений опреде-
т на каждом шаге безотносительно к значению на предыдущем.
Таким образом, ошибка не накапливается. Поскольку анализ жест-
копластического материала включает итерационную процедуру для
определения оптимального решения, продолжительность расчета
зависит от числа итераций на каждом шаге. Для уменьшения про-
должительности счета целесообразным оказывается метод, позво-
ляющий определить подходящую пробную скорость.
Контрольные вопросы
1. В чем суть метода конечных элементов для решения задач пластиче-
ского формоизменения?
2. Как производится построение дискретной модели для решения упруго-
пластических задач?
3. Как используется МКЭ для анализа жесткопластических деформаций?
133
Глава 6
МЕТОД ГРАНИЧНЫХ ЭЛЕМЕНТОВ
Метод граничных элементов нашел широкое применение для ре-
шения краевых задач во многих областях, в том числе механике
сплошной среды. Для любых краевых задач характерно наличие не-
которой области R, лежащей внутри границы С. Реальная задача
в области R моделируется дифференциальными уравнениями в част-
ных производных, решение которых отыскивается при определенных
ограничениях области. Если область R трехмерная, то С представля-
ет собой ограничивающую ее поверхность; в двухмерных задачах R —
плоская область, а С - ограничивающий ее контур.
В методах граничных элементов, в отличие от метода конечных
элементов, на элементы разбивается только граница С области.
Численное решение строится на основе полученных предваритель-
но аналитических решений для простых сингулярных задач
(например, задачи о сосредоточенной силе, приложенной в точке
бесконечной упругой среды) таким образом, чтобы удовлетворить
приближенно заданным граничным условиям на каждом элементе
контура С.
Поскольку каждое сингулярное решение удовлетворяет в R опре-
деляющим дифференциальным уравнениям в частных производных,
в этом случае нет необходимости делить область R на сетку элемен-
тов. В этом случае система уравнений, подлежащих решению, будет
значительно меньше системы метода конечных элементов.
Существует несколько методов граничных элементов: метод фик-
тивных нагрузок, метод разрывных решений, прямой метод гранич-
ных интегралов [5].
В настоящем учебном пособии, следуя [5], излагается только пря-
мой метод граничных интегралов для решения плоских задач.
6.1. Применение прямого метода граничных
интегралов для решения задач теории упругости
Рассмотрим плоскую краевую задачу теории упругости. Как из-
вестно, такие задачи характеризуются плоской областью R, ограни-
ченной контуром С. Область R может быть либо конечной (область
134
внутри контура С), либо бесконечной (область вне контура С), как
показано на рис. 6.1. В любом случае, с каждой точкой Q контура С
мы связываем касательные и нормальные смещения us и и и каса-
тельные и нормальные напряжения (или усилия) и ств. Эти величи-
ны задаются относительно локальной системы координат 5, п точки Q.
Если плоская задача теории упругости поставлена правильно, то
в каждой точке контура заданы касательное напряжение или каса-
тельное смещение и нормальное напряжение или нормальное смеще-
ние, т.е. две из четырех величин us, и, cts и ств должны быть заранее
известны. Оставшиеся две величины могут быть найдены как часть
решения задачи.
Рис. 6.1. Две области R для контура С: а - внутренняя область;
б - внешняя область
В этой главе для решения задачи мы используем прямой метод
граничных интегралов, который позволяет находить неизвестные
смещения и напряжения на границе прямо через заданные граничные
условия. В основе этого подхода лежит теорема линейной теории уп-
ругости, называемая теоремой взаимности.
6.2. Теорема взаимности и ее следствия
Теорема взаимности связывает решения двух различных краевых
задач для одной и той же области R. Теорема представляет прямое
следствие линейности уравнения равновесия и обобщенного закона
Гука.
Предположим, что одна краевая задача характеризуется смеще-
ниями us, и и напряжениями стд и ст на контуре С области R. Предпо-
135
ложим далее, что другая задача характеризуется смещениями
и', и'п и напряжениями а', с'„ на том же контуре С области R. То-
гда теорема взаимности утверждает, что работа, производимая пер-
вой системой сил (а5 и oj на перемещениях второй системы ( w’ и и' ),
равна работе, производимой второй системой сил (а' и <т'я ) на пере-
мещениях первой системы (и, и). Математически это можно запи-
сать в виде
)<*$ = J (<*'Л + .
(6.1)
где интегрирование выполняется по всему контуру С.
Если первая (нештрихованная) задача та, которую мы пытаемся
решить, а для второй (штрихованной) задачи мы уже знаем решение,
то (6.1) представляет собой интегральное уравнение, связывающее
неизвестные граничные параметры рассматриваемой задачи с задан-
ными граничными параметрами и с решением другой задачи для той
же области R. Решение этой второй задачи будем называть кон-
трольным, или тестовым, решением.
Предположим теперь, как и в ранее обсуждавшихся методах гра-
ничных элементов, что контур С можно аппроксимировать с помо-
щью N примыкающих друг к другу прямолинейных отрезков. Тогда
уравнение (6.1) можно представить в виде
JV N г
22 J №=22 J +ст
>=1 >=1 as'
(6-2)
где Д57 есть j-й отрезок, длина которого составляет 2а'. Если пред-
положить далее, что для рассматриваемой задачи смещения
и напряжения на границе в пределах каждого отрезка постоянны,
выражение (6.2) принимает вид
]<dS =
AS'
n'—’f
(63)
где aJs и cJn, и и',
отрезка.
7 - значения напряжений и смещений в центре J-ro
136
Поскольку число отрезков (граничных элементов) равно N, в ито-
ге мы имеем 4N граничных параметров и/, uJn, и . Половина
из них задана граничными условиями, которым мы пытаемся удовле-
творить, а другая половина - неизвестные, которые мы стараемся
найти. Предполагая пока, что нам также известно решение тестовой
задачи (u's, и'„, а'3, о'п), видим, что уравнение (6.3) содержит 2N не-
известных. Следовательно, для отыскания этих неизвестных необхо-
димо иметь еще 2N - 1 уравнений, аналогичных (6.3). Другими сло-
вами, в целом для заданной области R необходимо иметь 2N различ-
ных контрольных решений.
В следующем параграфе будет показано, что нужную систему 2N
контрольных решений можно найти, прикладывая к границе в цен-
тре каждого i-го элемента сосредоточенную касательную и нормаль-
ную силу F, и F„'. Считая, что N контрольных решений отвечают
сосредоточенным касательным силам F‘, i = 1, ..., N, запишем (6.3)
в виде
= Г <*', И ) dS J (f/ ) dS. (6.4)
as1 &sJ
Аналогично для N контрольных решений, отвечающих сосредо-
точенным нормальным силам F‘, i = I,...» N, имеем
Л
-,=' AS-*
Уравнения (6.4) и (6.5) можно представить в форме
(6.6)
137
где i принимает значения от 1 до N, а граничные коэффициенты
влияния (при обозначении коэффициентов влияния используются для
удобства такие же символы, как в методах фиктивных нагрузок
и разрывных смещений. Конечно, величины этих коэффициентов
в каждом случае различны)
dS,...
можно точно вычислить по контрольным решениям (см. п. 6.4).
Уравнения (6.6) можно преобразовать таким образом, чтобы 2N
неизвестных граничных параметров остались в одной стороне ра-
венств, а известные граничные параметры перешли в другую сторо-
ну. Если эти уравнения линейно независимы, они могут быть разре-
шены относительно неизвестных граничных параметров точно так
же, как решались системы уравнений в обсуждавшихся ранее методах
граничных элементов. Единственное отличие состоит в том, что (6.6)
позволяет нам найти неизвестные граничные смещения и напряжения
прямо по заданным граничным условиям.
6.3. Выбор контрольных решений
Представим бесконечную плоскость с контуром С (рис. 6.1). Если
в некоторой точке р приложить силу, то, используя результаты 4.2,
можно вычислить смещения и напряжения в любой точке плоскости.
В частности, можно вычислить касательные и нормальные смещения
и напряжения во всех точках контура С. В случае, когда область R
внутренняя, а точка р лежит вне контура С, из способа определения
напряжений следует, что если удалить область, внешнюю по отноше-
нию к С, то оставшаяся область, т.е. область R, будет находиться
в равновесии под действием напряжений, которые мы вычислили на С.
Следовательно, эти напряжения и соответствующие им смещения
образуют возможное для рассматриваемой задачи контрольное ре-
шение.
Это рассуждение не справедливо, если точка р лежит внутри кон-
тура С. Вводя в тело сосредоточенную силу, мы должны представить
вокруг точки малое отверстие и считать, что сила действует на кон-
туре этого отверстия. Тогда область R имеет две границы - контур С
и границу отверстия вокруг точки р. Поэтому, если рассматривать
138
только напряжения, приложенные на границе С, они не будут урав-
новешены и не дадут приемлемое контрольное решение.
Проведенные рассуждения остаются в силе и в случае, когда рас-
сматриваемая область R представляет неограниченную область вне
контура С. Следовательно, для обеих задач - и внешней и внутренней
(рис. 6.1) - необходимо, чтобы при отыскании контрольного решения
сосредоточенная сила была приложена в точке р вне области R.
Фактически сосредоточенная сила обеспечивает два контрольных
решения, поскольку она разлагается на две составляющие, образую-
щие между собой прямой угол. Каждая составляющая обеспечивает
решение для смещений и*,и'п, напряжений <^, сг„, которые, будучи
проинтегрированы согласно (6.3), дают коэффициенты в системе (6.4)
и (6.5). Для получения коэффициентов 2W уравнений мы вводим N
сил в N различных точках неограниченной плоскости, но не в самой
области Л. За исключением этого единственного ограничения, точки N
могут быть выбраны произвольно, но при этом нет гарантии, что
полученная система уравнений будет хорошо обусловленной или да-
же линейно независимой. Последовательный подход, который, как
установлено, приводит к хорошо обусловленным уравнениям, за-
ключается в выборе N контрольных точек в серединах W отрезков
контура С. Точнее, мы прикладываем сосредоточенные силы к сере-
динам N отрезков снаружи области R.
6.4. Коэффициенты влияния
В прямом методе граничных интегралов граничные коэффициен-
ты влияния получают путем приложения сосредоточенной силы
(с компонентами Fj и F’ ) в средней точке г-го отрезка контура С
и интегрирования смещений и напряжений, вызванных этими сила-
ми, вдоль у-го отрезка в соответствии с (6.4) и (6.5). Изменяя i от 1 до
N, т.е. прикладывая поочередно N сосредоточенных сил по контуру,
получаем необходимую систему алгебраических уравнений (6.6).
При вычислении коэффициентов влияния В‘£, ... в (6.6) для удоб-
ства воспользуемся локальной системой координат х, у с началом
в центреу-го отрезка контура С (рис. 6.2). Эти координаты соответст-
вуют обычным локальным координатам $ , пу, причем ось х (или s')
совпадает с направлением обхода контура, а ось у (или п j направ-
лена вне рассматриваемой области (рис. 6.1). Мы хотим вычислить
139
смещения и напряжения на у-м отрезке, вызванные действием сосре-
доточенной силы, приложенной в центре i-ro отрезка (в точке [j] на
рис. 6.2). Компоненты этой силы, параллельная и перпендикулярная
i-му отрезку, равны F* и F' соответственно. Из рис. 6.2 легко видеть,
что компоненты силы в направлениях х = sJ и у = nJ определяются
следующими выражениями:
= Ff'cos у - F„ sin у,
F-y = F‘ sin у + F„' cos y,
(6.7)
где y = p/-p/.
Используя результаты 4.2, можно записать выражения для смеще-
ний и напряжений, вызванных силой F*, Fy, приложенной в точке
х = = Су,т.е. в точке [i] (рис. 6.2). Это осуществляется просто пу-
тем замены х и у в (6.7)-(6.8) на х = сг, у - с у и, конечно, изменени-
Рис. 6.2. Определение локальных координат для граничных коэффициентов
Тогда смещения и'3 = и*,и'п = и- и напряжения а ’3 = <3уу,а’п =Оуу
на у-м отрезке (наше контрольное решение) находим из полученных
выражений, подставляя в них J = 0. Используя величины Fj,Fy из
(6.7), находим смещения для контрольного решения
140
J
2GL
2GL
9 — s
и
2G
(3-4v)gsiny-
—c^cosy
dx
dg
(6.8)
dg
* n
2G
и напряжения
где функция g(x, j) имеет вцд:
5
g(x,y)=----т—-ln[(x-c7)2 +(.и-с7)2]2. (6.10)
4*(1- v)
Сама функция и ее производные определяются в (6.8) и (6.9) при
J = 0.
Как отмечалось ранее, выражения (6.8) и (6.9) дают фактически
два контрольных решения для каждого элемента i: одно для каса-
тельной силы F' , а другое для нормальной силы F' . Коэффициен-
ты влияния в (6.6) получаются путем поочередного выбора
141
этих решений (т.е. F' О, F' = О и F* — О, F* * 0), подстановки (6.8)
и (6.9) в (6.4) и (6.5) и выполнения необходимого интегрирования
вдоль As7, т.е. интегрирования по х в пределах от - а7 до + а7. Вы-
полнив указанные операции, видим, что нет нужды в определении
величин F* и F *, поскольку они входят как сомножители в обе части
уравнений. Следовательно, можно считать, что F' и F' равны +1.
Поступая таким образом, приходим к следующим результатам:
+aJ / 1
Bfs= I ---[(3-4v)7}cosy + cf(t^ sin у-T3 cos у)],
Jj 2G
M = f и'я(л,' )dx =------[(3-4v)7'i sinY + OyCr? cosy+ ?з siny)],
J, 2G
“fl7
+fly J
= f и'АВп)^х=--------[-(3-4v)t; siny + Cj; (T2 cosy + T, siny)}
J, 2G
-a*
2G
[(3-4v)Tj cosy-Cj; (Tj siny-T^ cosy)]; (6.11)
+c? i _ _ _ _
I n's(F/)dx =—[(l-2v)T2 siny+2(1-v)T3 cosy+с5(т4 siny+Г5 cosy)}
2G
—aJ
л = [-0-2v)T2cosy+2(l-v)7’3 siny+c- (t4 cosy-T5 siny)}
= jo'i(F^)aK=[(l-2v)7;2cosy-2(l-v)T3siny+Cp(7;cosy-r5siny)],
+a^
= |o,„(F^)c^=[(l-2v)T2siny+2(l-v)T3cosy-c5i(T4siny+T5COsy)]. (6.12)
—fl-'
В этих формулах величины 7],..., Т5 представляют собой опреде-
ленные интегралы от функции g(x, у) и ее производных
(вычисленных при у — 0 ). Они равны
142
Выше при обсуждении мы неявно предполагали, что точки [ij и [/]
на рис. 6.2 не совпадают, т.е. величины с учетом (6.13) прием-
лемы и в случае совпадения точек, если при вычислении многознач-
ной функции арктангенса использовать обычную процедуру пре-
дельного перехода. Подставляя с- = 0 и переходя к пределу при с?,
стремящемся к нулю со стороны положительных значений у (т.е. из-
вне области R, как описано в (6.3), находим, что в (6.13) отличны от нуля
только величины Т], Т'з . Эти величины определяются выражениями
(6-14)
Диагональные члены коэффициентов влияния, отражающие соб-
ственные воздействия элементов, получаются использованием этих
результатов с (6.11) и (6.12) с учетом того, что у = р' -0^ — 0:
143
^=^=0, В“=В''п =---------3 4v у Ina';
4nG(l-v)
(6.15)
(6.16)
6.5. Построение системы уравнений
Уравнения (6.6) составляют основу прямого метода граничных
интегралов. Для любой краевой задачи половина из 4У параметров
этих уравнений (т.е. для i = 1, ..., N) задаются как гра-
ничные условия, в то время как другая половина соответствует неиз-
вестным. Коэффициенты влияния определяются в соответствии
с геометрией задачи по формулам (6.11) и (6.12). Следовательно,
уравнения (6.6) можно использовать для записи системы 2N алгеб-
раических уравнений с 2N неизвестными точно так же, как это дела-
лось в рассмотренных выше методах граничных элементов. Неиз-
вестными в этой системе уравнений являются фактические граничные
смещения или напряжения, которые не заданы как граничные условия.
В краевой задаче в напряжениях, например для всех N элементов,
известны величины ст' = (о^)0 и <з‘п = (ctJJ0. Тогда левая часть урав-
нений
(6.17)
оказывается известной (ср. (6.6)), а смещения u‘s,и‘п, / -1,..., N, можно
найти обычным способом, решая эти уравнения. Аналогичную сис-
тему уравнений можно построить и в случае, когда на всех N элемен-
тах известны смещения м* = )о и ип = (wn )о» тогДа роль неизвест-
ных играют напряжения ст, и ctJ, (/ = 1, ...,А).
Для смешанной краевой задачи система уравнений строится путем
преобразования (6.6) таким образом, чтобы 2А известных граничных
параметров были в одной части равенств, a 2N неизвестных пара-
метров - в другой. Например, если на всех N граничных элементах
144
заданы нормальные напряжения <т'„ = jj и касательные смещения
«' = (и' )0, систему уравнений можно записать в виде
(6.18)
В этой системе неизвестными теперь служат касательные напря-
жения Oj и нормальные смещения и'п для i = 1,...» N. Следует также
заметить, что запись (6.18) неудобна, поскольку коэффициенты влия-
ния (В^ и В£) содержат множитель (2G)~l, и, следовательно, они
обычно на несколько порядков меньше, чем коэффициенты смеще-
ний (А% и (ср. (6.11) и (6.12)). Если коэффициенты при напря-
жениях умножить на постоянную 2G, а сами напряжения разделить
на ту же постоянную, уравнение (6.18) можно переписать в виде
(6-19)
где теперь все коэффициенты влияния имеют одинаковый порядок.
Аналогичный способ можно использовать для других типов смешан-
ных краевых задач, включая те, в которых на одной части границы
заданы напряжения, а на другой - смещения.
Согласно предыдущим рассуждениям, основную систему алгеб-
раических уравнений для всех типов краевых задач можно предста-
вить в виде
(6.20)
145
В этих уравнениях У/ и У„, 1=1,N, суть определенные линей-
ные комбинации известных граничных параметров, а ... - коэф-
фициенты влияния при неизвестных граничных параметрах
X/ и XJn . Например, коэффициент равен , если на J-м эле-
менте неизвестно касательное вмещение и/ = но он равен
- (2<7)в£, если на этом элементе неизвестно касательное напряжение
(2Gr’a;(r/).
6.6. Формулы Сомильяны
До сих пор прямой метод граничных интегралов мы привлекали
только для вычисления неизвестных смещений или усилий на границе
С произвольной области R. Если же мы хотим найти решение внутри
рассматриваемой области R, можно воспользоваться интегральными
тождествами, известными как формулы Сомильяны. Для плоской
задачи эти формулы дают смещения внутренней точки р области R
в виде
«х(Р) = -J +
+J \ )+«»«»(Л:
+ (Л )+ (Fy )]<*£,
с
(6.21)
где и, ия и а, сп - граничные смещения и напряжения для задачи,
решение которой мы уже имеем, a u's (Ft), и’п (F() и o', (F, ), ст^, (Fz ) -
нормальные и касательные смещения и напряжения на границе С,
вызванные действием сосредоточенной силы F,- =(FX, Fy )=(+!,+ !)
в точке р. У читателя может возникнуть искушение приравнять пра-
вые части (6.21) нулю на основании теоремы взаимности. Этого
нельзя делать, поскольку теперь мы приложили сосредоточенную
силу в точке, лежащей внутри области R. Поэтому из тех рассужде-
146
ний, которые приводились в 6.3, следует, что теорема взаимности
применима не к границе С, а к границе, которая помимо С включает
малый контур вокруг точки р. Формулы Сомильяны действительно
получают, используя теорему взаимности для такой границы
и стягивая точку в малый контур вокруг р.
Уравнения (6.21) можно разрешить численно, разбивая границу С
на элементов и предполагая, как и ранее, что us, ип и ат, оп постоян-
ны в пределах каждого элемента. Это дает
(6.22)
где теперь все граничные параметры us,un hos, известны. Инте-
гралы в (6.22) являются коэффициентами, которые необходимо вы-
числить для каждой точки р, где отыскиваются значения смещений
«х(р) и и (р). Дифференцируя выражение для смещений в этих точках,
находим деформации. Далее, используя соответствующую форму
обобщенного закона Гука, определяем напряжения в точке р. Выра-
жения для «внеграничных» коэффициентов влияния смещений и на-
пряжений даются ниже.
6.6.1. Коэффициенты смещений
Интегралы в (6.22) вычисляются в локальной системе координат
х,у с началом в центре у-го элемента рис. (6.3).
Точка р на рисунке является произвольной точкой с координата-
ми х - с*, у = Су, не лежащей на границе С.
147
Рис. 6.3. Определение локальных координат для коэффициентов
в точках вне границы
Штрихованные величины в (6.22) представляют смещения
u's =их, и'п —Uy и напряжения ст, = ст , <з'„ — стуу, возникающие в J-m
граничном элементе под действием сосредоточенных сил F — +1,
F = +1, приложенных в точке р. Выражения для этих величин можно
получить непосредственно из (6.8) и (6.9), полагая у = -р-' (т.е. pJ = 0)
и F* = Fx, F„ = Fy. В (6.22) необходимые интегралы находятся точно
так же, как ранее, путем рассмотрения в отдельности случаев F' = +1,
F = +1. Уравнения (6.22) тогда можно записать в виде
^[-(3-4v)Tjsinpy
2G
4g
(6.23)
148
где Г) и другие коэффициенты определяются по (6.13).
6.6.2. Коэффициенты напряжений
Формулы для напряжений в точке р можно найти, подсчитав де-
формации, отвечающие смещениям (6.23), и используя далее обоб-
щенный закон Гука для случая плоской деформации. Деформации
даются соотношениями:
ОС
* ^yy(P) д ’ ^Ху^~РУ д
ОС,, 2
(6-24)
где сх и су - компоненты х и у вектора, соединяющего центр /-го гра-
ничного элемента с точкой р (рис. 6.3). Используя формулы преобра-
зования
C-s = с„ cos
0'
+ casing7;
с-у = -сх sin Р7 + су cos р7,
(6.25)
находим по цепному правилу частного дифференцирования
д 8с- 8 де? д j д . ni д
—— = —Л-----+ — -----= cosp7-----smB7-----;
дсх дсх 8сх 8сх дсу дсх 8су
8 8сх 8 tdc~y 8 . oi 8 of 8
---= —------+—— = smB7 +cosp7-.
8cv 8cv 8c- 8cv дс~-----------------------дс?-дс7
Следовательно, деформации (6.24) можно записать в виде
(6.26)
149
, . диАр) „j диАр) . .
ехх(р) = —~ cosp7-------р7;
дс* дст
Suv(p) , диАр)
еуу{р) = —г----sin р' + —~ cos ру;
осу дс„
-У
, ч 1 ди Ар) . л. диАр) ni
еху(Р) = 7 —Т-----Sin р7 + —-cos р7 +
2 дс^ дсъ
(6.27)
диу{р) aJ ди (р) . .
—-----cos В7--------smB7
дс— дс^
л л
Подставляя (6.23) в (6.27), находим искомые формулы для напря-
жений. Окончательные выражения имеют вид
XX
5 sin2pJ -
N г _ , _ \1
+ 2<?^Г75 -сДГ^тгр-' +T7cos2p7 )]u' +£[-7’2 -2(1-v)x
/- - • (- -М
х^Т2 cos2p7 -7*з sin2p7 )+Су (Г4 cos2p7 +7*5 sin2p7 J] а7 +
+ £ [- Hl ~ v )(r2 sin 2p7 + 7*3 cos 2p7 )]+cy (t4 sin 2p7 - T5 cos2 p7 )] o7 ,
7=1
'7 sin2p? )]us
'6 sin2pJ + T7 cos2py J
UJ
x (7*2 cos2p7 - T3 sin2p7
(6.28)
s **"
n >
150
В этих выражениях функции 7],Т5 даются формулами (6.13),
а функции Т6, Tj определяются выражениями
6.7. Преобразование координат
До сих пор в данной книге все результаты выражались через ло-
кальные касательные и нормальные компоненты смещений и усилий
вдоль заданной границы. Однако в большей части литературы мето-
ды граничных элементов формулируются в терминах глобальных
компонентов х и у этих величин, т.е. u. - (ux, и t. = (t, t}. Связь
между локальной и глобальной формулировками метода граничных
элементов легко установить, используя простые формулы преобразо-
вания координат. Например, используя соотношения
и{ = ul cosp7 + id sin в7,
uJn = -tdx sin p7 + uj, cosp7;
cd = d cos p7 + d sin p7,
•3 A J *
o7 = -d sin P7 + d cos p7,
(6.30)
(6.31)
уравнения (6.6) можно представить в виде
If N N
(6-32)
151
где «глобальные» коэффициенты влияния В%х... связаны с локаль-
ными коэффициентами выражениями вида
В”х = Bl cosру - В*п sinрЛ .... (6.33)
Необходимые результаты можно также получить, не обращаясь
к локальным координатам. Проиллюстрируем это ниже для прямого
метода граничных интегралов. С целью ознакомления читателя
с терминологией, используемой в литературе, часть вывода проведем
в индексных обозначениях.
Начнем с записи теоремы взаимности (6.1) в виде
+ и/у ~}dS.
(6.34)
Здесь, как и прежде, штрихованные величины относятся к нашему
контрольному решению. В индексных обозначениях это соотношение
имеет вид
[ t,u' dS = j dS, (6.35)
где по повторяющемуся координатному индексу производится сум-
мирование.
Следуя далее [6], рассмотрим две произвольные точки на конту-
ре С. Одну из этих точек назовем «точкой поля», а другую - «точкой
нагружения». Наше контрольное решение отыскивается путем рас-
смотрения смещений и усилий в точке поля, вызванных действием
сосредоточенной силы в точке нагружения. (Предположим пока, что
эти точки не совпадают). Как и ранее, заметим, что сосредоточенная
сила в точке нагружения порождает два контрольных решения, по
одному для каждой компоненты силы. Это позволяет получить вы-
ражения для обоих решений одновременно.
Если точка нагружения имеет относительно точки поля коорди-
наты х — с , то смещения для контрольного решения можно записать
в виде
152
где g(Xj, *2) дается формулой
(6.37)
Используя индексные обозначения, уравнения (6.36) можно запи-
сать в компактной форме:
I
где U - «тензорное поле», определяемое по формуле
(6.39)
2G
Усилия t для контрольных решений получаем из соотношений
(6.40)
где n'j задает компоненты единичной внешней нормали к границе С в точ-
ке поля. Напряжения а'7 = , в свою очередь, находим, используя
соотношения
У
где деформации e'tJ связаны со смещениями (6.38) определениями
е'у = — +ujj). Читателю следует убедиться, что результаты
вычисления можно представить в виде
ЭТОГО
(6.41)
где
8 +
О Дг + —----
Подстановка (6.41) в (6.40) дает:
I
(6.43)
153
или после суммирования по повторяющемуся индексу
4 ~ (^1ЛИ1 + $2ikn2 )^к •
(6.44)
В этом выражении повторяющийся индекс к действует как немой
и может быть заменен любым символом без изменения смысла выра-
жения. В частности, (6.44) можно записать в виде
где Т - «тензорное поле»:
(6-45)
(6.46)
Если точка Q границы С есть точка поля, а Р - точка нагружения,
то (6.38) и (6.45) можно записать следующим образом:
и! е = ^(Л(2)Г/Р);
/;с=т7;(ле)^(Р).
(6-47)
Рассматривая отдельно два контрольных решения, представлен-
ных условиями
FJP) * О, Г2(Р) = 0 и Fj(P) = О, F2(P) * О,
воспользуемся выражениями (6.47) и соотношением взаимности (6.35)
для получения следующих двух уравнений:
Z U
р.ЧСХЛЛЛС)^ Р)Л(б) =(2)^ (Р,е)Д(РЖб);
г с. (6.48)
г =[ н,(№2/(ле)^2(Р)Л(е),
или, поскольку F((P) и F2(P) произвольны,
J ti(.QWP,Q)ds(Q) =| «ЛОТиСЛСже);
J =f «.(№(ле>ад).
(649)
Заметим, что оба этих соотношения можно получить из равенства
J tt(QWfi(P,Q)d^Q) =f ц(2)7}ЛЛе)Л(СХ (6.50)
154
принимая j = 1 и J = 2. В случае,когда точка поля Q и точка нагруже-
ния Р совпадают, необходимо рассматривать предельный переход,
в котором сосредоточенные силы FX(P) и F2(P) приложены к границе
извне рассматриваемой области R (п. 6.3).
6.7.1. Численное решение
Соотношения (6.49) можно записать в развернутой форме:
|к(е)^х(ле)+^(е)^(лок(0)=
= jk (Q)TXX (Р ,Q) + иу (Q)Txy (P,Q)]ds(Q\,
j[rx (Q)Uyx (P ,Q) + ty (Q)Uyy ( P,e)>(6) =
= Jk (2)7,, (Л6) + uy (Q)Tyy (P.Q)] ds(Q).
(6.51)
При численном решении этих уравнений граница области С, как
и ранее, аппроксимируется N прямолинейными отрезками, примы-
кающими друг к другу, и предполагается, что смещения и усилия
в пределах каждого отрезка постоянны. Поступая так же, как
и в п. 6.2, получаем
В этих уравнениях верхние индексы i и j относятся к номерам гра-
ничных элементов, и их не следует путать с аналогичными символа-
ми, использованными ранее в качестве координатных индексов. Если
не обращать внимания на очевидные различия в обозначениях, то
(6.52) совпадает с (6.32).
155
6.7.2. Формулы Сомильяны
В заключение этого раздела заметим, что формулы Сомильяны
(6.21) также можно записать через компоненты х и у граничных сме-
щений и усилий. Точка нагружения р в этом случае лежит внутри
рассматриваемой области. Используя обозначения, принятые ранее,
имеем
«!</>)=-J «,(e)75/(Ae)^(e)+p,(eWu(Aewe);
с с
«2 (Р) = Ч;(в)Г2.( Р, Q)ds{Q) +| /, (e)l/2J( p,Q)ds(Q),
с с
(6.53)
или кратко:
«//>)=-J ц (е)гу/( ле)<*«2) +J '«(ехшоад). (6.54)
где i и j принимают значения 1 или 2.
Расчеты по (6.53) можно выполнить численно, разбивая границу
С на N элементов и считая, что смещения и и усилия в пределах ка-
ждого элемента постоянны. Тогда функции Т (р, Q) и U (p, Q) можно
проинтегрировать для элементов, определяя тем самым граничные
коэффициенты влияния для смещений в точке р. Напряжения в этой
точке определяются с помощью процедуры, описанной в (6.7).
6.8. Применение прямого метода граничных
интегралов для решения задач вязкоупругости
6.8.1. Вязкоупругое поведение материала
Упругие тела и вязкие жидкости существенно различаются своими
свойствами при деформировании. Упругие деформируемые тела после
снятия приложенных нагрузок возвращаются к своему естественному,
или недеформированному, состоянию. В отличие от них несжимаемые
вязкие жидкости совсем не имеют тенденции возвращаться после сня-
тия нагрузки в исходное состояние. Кроме того, напряжения в упру-
гом теле связаны непосредственно с деформациями, в то время как
напряжения в вязкой жидкости зависят (за исключением гидростати-
ческой составляющей) от скоростей деформации.
Поведение материала, которое объединяет в себе оба этих свойст-
ва - и упругости, и вязкости, - называют вязкоупругим. Упругое тело
156
и вязкая жидкость занимают крайние противоположные точки в ши-
роком спектре вязкоупругих сред. Хотя вязкоупругие материалы чув-
ствительны к температуре, последующее изложение ограничивается
условиями изотермин и температура входит в уравнения только как
параметр.
6.8.2. Простейшие механические модели
вязкоупругого поведения
Линейную вязкоупругость для одномерного состояния удобно
трактовать при помощи механических моделей, которые наглядно
демонстрируют поведение различных вязкоупругих материалов. Эти
модели строятся из таких механических элементов, как линейно-
упругая пружина с модулем упругости G (массой этой пружины пре-
небрегают) и вязкий элемент (демпфер) с коэффициентом вязкости
(вязкий элемент представляет собой поршень, движущийся в цилинд-
ре с вязкой жидкостью). Как показано на рис. 6.4, сила ст, растяги-
вающая пружину, связана с ее удлинением а формулой
ст = Ge . (6.55)
Подобное же соотношение существует и для демпфера
о = Т]€ , (6.56)
где ё = de/dt. Можно придать большую общность этим моделям
и устранить размерные эффекты, если в качестве ст рассматривать
напряжение, а в качестве е - относительную деформацию.
Рис. 6.4. Механические модели:
а - линейный упругий элемент; б - вязкий элемент
157
Модель Максвелла вязкоупругого тела является комбинацией
пружины и вязкого элемента (демпфера), соединенных последова-
тельно (рис. 6.5, а). Модель Кельвина или Фойхта представляет собой
параллельное соединение тех же элементов (рис. 6.5, б).
Соотношение между напряжением и деформацией (фактически
содержащее также и их скорости) для модели Максвелла дается фор-
мулой
- + - = £, (6.57)
G т|
а для модели Кельвина формулой
а = Ge + г|ё.
а
Рис. 6.5. Механические модели: а - модель Максвелла; б - модель Кельвина
(6.58)
Эти уравнения являются по существу определяющими уравне-
ниями вязкоупругости в одномерном случае. Полезно написать их
в операторной форме, используя линейный оператор дифференцирова-
ния по времени д = dfdt. Тогда уравнение (6.57) будет иметь вид
(6.59)
6.8.3. Трехмерная теория
При создании трехмерной теории линейной вязкоупругости
обычно принято рассматривать отдельно вязкоупругое поведение
в условиях так называемого чистого сдвига и чистого расширения [6].
Таким образом, эффекты искажения формы и изменения величины
объема изучаются независимо и затем их описания комбинируются,
чтобы построить общую теорию. Математически это обеспечивается
разложением тензоров напряжений и деформаций на их девиаторную
и шаровую части, для каждой из которых затем пишутся определяю-
158
щие соотношения вязкоупругости. Разложение тензора напряжений
рмулой
дано
^ij^kk >
(6.60)
а разложение тензора малых деформаций - формулой
&У ~^^у^кк /^'
(6.61)
Принимая обозначения, использованные в этих равенствах, запи-
шем в операторной форме трехмерное обобщение определяющих
уравнений вязкоупругости (6.64):
{Р}^ = 2 {Q} е(. (6.62а)
и
{Л0а/у = 3{Ме„, (6.626)
где {Р}, {Q}, {Л/}, {TV} - дифференциальные операторы, определен-
ные в соответствии с (6.65), а числовые множители введены для удоб-
ства. Так как практически все материалы упруго реагируют на уме-
ренные гидростатические нагрузки, в качестве операторов {М}
и {7V}, связанных с расширением, обычно берут постоянные и преоб-
разуют соотношения (6.62) к виду
{РЦ. = 2{2}е?, (6.63а)
ай=ЗХе„, (6.636)
где К - объемный модуль упругости.
Линейное дифференциальное уравнение можно записать символи-
чески
{PJo = {Q}e,
(6.64)
где операторы {Р} и {Q} определены формулами
Q = G + r\dt.
(6.65)
Для модели Кельвина
Р = 1;
(6.66)
6.8.4. Анализ вязкоупругого напряженного состояния
Принцип
тветствня. Задача исследования напряженного со-
стояния в некотором изотропном вязкоупругом теле, которое зани-
159
мает объем К и ограничено поверхностью S (рис. 6.7), ставится сле-
дующим образом. Пусть заданы массовые силы Ь., действующие
внутри V, и пусть на части S, границы тела известны внешние по-
верхностные силы tfn)(xk,l), а на части S2 поверхности тела - смеще-
ния поверхности gt(xk, t).
Рис. 6.7. Схема нагружения изотропного вязкоупругого тела
Тогда система, дающая постановку задачи, включает следующие
соотношения:
1) уравнения движения (или равновесия)
а,7у+6,- = р«,;
(6.67)
2) соотношения, выражающие деформации через перемещения
2 6^ = (л,у +«;,), (6.68)
или скорости деформации через скорости
2ё,7 =(иг7+Оу,); (6,69)
3) граничные условия
,t) на Sj, (6.70)
на S2> (6-71)
160
4) начальные условия
м,(хк,О) = мо,
wffa>0) = vo;
5) определяющие уравнения, записанные
(6-72)
(6.73)
через линейные ди
ренциальные операторы, т.е. в форме (6.65), (6.66).
Если форма тела и условия нагружения достаточно просты и если
поведение материала может быть представлено одной из простейших
моделей, то приведенную выше систему уравнений можно проинтег-
рировать непосредственно. Однако для более общих условий обычно
принято искать решение, пользуясь принципом соответствия упругой
и вязкоупругой задач. Этот принцип основывается на том, что система
основных уравнений теории упругости и преобразования Лапласа по
времени вышеприведенной системы основных уравнений теории вяз-
коупругости записываются одинаково. Соответствующие уравнения
для квазистатических изотермических задач, в которых черточки оз-
начают преобразования Лапласа по времени, например,
ОО
f{xk’s
О
(6.74)
сопоставлены в табл. 6.8.
Таблица 6.8
Сопоставление основных уравнений теории упругости и вязкоупругости
Уравнения теории упругости
ст..у+/,. =0
2e,j — \Uy )
<3tjnj=t\n^ на 5,
Uf = g, на S2
Sy ~ 2Gey
C„ = 3K£lt
Преобразования уравнений теории
вязкоупругости
°ij,j'+bi =0
Gytij на
й( = g, на S2
P(s)Sij = 2Qeij
Og — 3-КЁц
Из этой таблицы видно, что если G в уравнениях теории упруго-
сти заменить на QIР, то обе группы уравнений будут иметь одина-
ковую форму. В силу этого, если в решении соответствующей задачи
теории упругости G заменить на Q! Р для вязкроупругого материа-
ла, то полученный результат будет преобразованием Лапласа реше-
161
ния задачи теории вязкоупругости. Обратным преобразованием най-
дем само решение для вязкоупругого материала.
Этот принцип соответствия может быть установлен и для задач,
отличных от квазистатических.
Контрольные вопросы
1. В чем суть метода граничных элементов?
2. В чем различия в МКЭ и МГЭ?
3. В чем особенности применения МГЭ для решения задач вязкоупру-
гости?
Заключение
Данное учебное пособие позволит будущим инженерам освоить
методы математического моделирования процессов обработки ме-
таллов давлением, поставить соответствующую краевую задачу,
сформулировать начальные и граничные условия, выбрать наиболее
эффективный метод решения.
Решение краевой задачи дает возможность определять в области
пластического течения распределение скоростей, деформаций, темпе-
ратур, напряжений, использование ресурса пластичности.
В связи с этим появляется возможность анализа условий
заполнения профиля в процессах пластического формоизменения.
Это позволяет находить предельную степень деформации, работу
деформации, местный перегрев металла, вызванный локальным
выделением тепла, определять нагрузки на инструмент, оценивать
его износ, прочность, количество переходов, ресурс пластичности,
прогнозировать механические характеристики и свойства изго-
тавливаемых деталей.
Появляется возможность оптимизации процесса с точки зрения
энергосиловых затрат, стойкости инструмента, расхода металла, ка-
чества получаемых изделий.
Эти знания необходимы для анализа существующих процессов
обработки металлов давлением с целью их совершенствования и оп-
тимизации, а также для разработки новых высокоэффективных ре-
сурсосберегающих технологий.
163
Список литературы
1 . Качанов JLM. Основы теории пластичности. - М.: Наука, 1969. -420 с.
2. Гун ГЯ. Теоретические основы обработки металлов давлением. - М.: Ме-
таллургия, 1980. - 456 с.
3. Применение ЭВМ для решения задач теплообмена / Т.Н. Дульнев,
ВТ. Парфенов, А.В. Сигалов. - М.: Высш, ппс, 1990. - 207 с.
4. Теория пластических деформаций металлов / Е.П. Уиксов, У. Джонсон,
В.Л. Колмогоров и др. — М.: Машиностроение, 1983. -598 с.
5. Крауч С, Стпарфилд А. Методы граничных элементов в механике твер-
дого тела. - М.: Мир, 1987. - 328 с.
6. Дж. Мейз. Теория и задачи механики сплошных сред. ~ М.: Мир,
1974. - 320 с.
ОГЛАВЛЕНИЕ
Введение.......................................................3
Глава 1. Основные положения механики сплошных сред.............4
1Л. Напряженное состояние деформируемой среды..................4
1.2, Деформация...............................................11
13. Скорость деформации......................................16
1.4. Дифференциальные уравнения равновесия.
Граничные и начальные условия.................................21
1,5. Условия текучести...................................... 23
1,6. Соотношения между напряжениями, деформациями
и скоростями деформации.......................................24
1.7. Некоторые обобщения.................................... 34
1.8. Уравнение теплопроводности. Начальные
и граничные условия...........................................35
Глава 2. Теория разрушения и пластичность металлов............46
2,1. Гипотеза о разрушении металлов при пластической
дейювмации....................................................46
Глава 3. Математическое моделирование процессов обработки
металлов давлением.............................................52
3.1. Краевая задача обработки металлов давлением...............52
3.2. Методы решения краевых задач..............................59
Глава 4. Метод конечных элементов для решения задач
теплопроводности...............................................61
4.1. Основные концепции метода конечных
элементов (МКЭ).................................................61
4,2. Построение дискретной модели и функций
формы элементов.................................................65
4.3. Система уравнений метода конечных элементов.
Локальная и глобальная матрицы..................................69
4.4. Формирование глобальных матрицы и вектор-
столбца. Решение системы уравнений МКЭ..........................76
4.5. Программная реализация МКЭ................................82
Глава 5. Метод конечных элементов для решения задач
пластического формоизменения...................................92
5.1. Основные характеристики элементов.........................92
5.2. Использование МКЭ для анализа упругопластической
плоской деформации............................................ 102
5.3. Теория МКЭ для упругопластической деформации..............110
5.4. Использование МКЭ для анализа жесткопластических
деформаций................................................... 119
5.5. Применение МКЭ в процессах штамповки......................131
165
Глава 6. Метод граничных элементов..........................134
6.1. Применение прямого метода граничных интегралов для
решения задач теории упругости..............................134
6.2. Теорема взаимности и ее следствия.....................135
6.3. Выбор контрольных решении......................... 138
6.4. Коэффициенты влияния.............................. 139
6.5. Построение системы уравнений..........................144
6.6. Формулы Сомильяны................................... 146
6.7. Преобразование координат............................ 151
6.8. Применение прямого метода граничных интегралов
для решения задач вязкоупругости............................156
Заключение.................................................163
Список литературы..........................................164
Учебное издание
Илья Борисович Покрас
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
ПРОЦЕССОВ ОБРАБОТКИ МЕТАЛЛОВ ДАВЛЕНИЕМ
Учебное пособие
Гл. редактор (директор) Г.А. Осипова
Редактор Л/.А. Ложкина
Технический редактор Н.К. Швиндт
Корректор Н.К. Швиндт
Верстка Т.Н. Пермяковой
Подписано в печать 06.06.02. Формат 60x84/16. Бумага офсетная
Гарнитура Times. Усл. печ. л. 9,76. Уч.-изд. 7,92 л.
Тираж 200 экз. Заказ № 271
Отпечатано в типографии Издательства ИжГТУ
Издательство Ижевского государственного технического
университета. 426069, Ижевск, Студенческая, 7
И.Б. Покрас
МАТЕМАТИЧЕСКОЕ
МОДЕЛИРОВАНИЕ
ПРОЦЕССОВ ОБРАБОТК
МЕТАЛЛОВ ДАВЛЕНИЕМ