Вариационная постановка и разработка методов решения задач контактного взаимодействия тел при конечных деформациях тема автореферата и диссертации по механике, 01.02.04 ВАК РФ

UUJi 72298

На правах рукописи

МОРЕВ ПАВЕЛ ГЕННАДЬЕВИЧ

ВАРИАЦИОННАЯ ПОСТАНОВКА И РАЗРАБОТКА МЕТОДОВ РЕШЕНИЯ ЗАДАЧ КОНТАКТНОГО ВЗАИМОДЕЙСТВИЯ ТЕЛ ПРИ КОНЕЧНЫХ ДЕФОРМАЦИЯХ

Специальность 01 02 04 - Механика деформируемого твердого тела

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

1 g

Тула 2008

003172298

Работа выполнена в Орловском государственном техническом университете

Научный руководитель доктор физико-математических наук,

профессор

Маркин Алексей Александрович

Официальные оппоненты доктор физико-математических наук,

профессор

Бровко Георгий Леонидович

доктор физико-математических паук, профессор

Шоркин Владимир Сергеевич

Ведущая организация Пермский Государственный технический

университет

Защита диссертации состоится "01" июля 2008 г в 10 — часов на заседании диссертационного совета Д 212 271 02 при ГОУ ВПО "Тульский государственный университет" по адресу 300600, г Тула, пр Ленина 92 (9101)

С диссертацией можно ознакомиться в библиотеке ГОУ ВПО "Тульский государственный университет"

Автореферат разослан " 70 " мая 2008 г

Ученый секретарь -

диссертационного совета JIА Толоконников

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы Сейчас трудно представить создание новых производственных технологий без разработки адекватных математических моделей и их всестороннего исследования в численных экспериментах Для технологий, связанных с формоизменением твердых тел, конечным этапом этих теоретических исследований является решение соответствующей краевой задачи, которая, как правило, оказывается контактной Вместе с тем такие задачи во многих отношениях - одни из самых сложных краевых задач математической физики Ввиду своей важности и сложности контактные задачи привлекали большое число исследователей как в нашей стране (А С Кравчук, В С Давыдов, Е Н Чумаченко, Э Р Гольник, Н И Гундорова, А А Успехов и др ), так и за рубежом (G Pietrzak, A Curnier, F Armero, Е Petoch, Р Alart, М Barboteu, F Lebon, D Barlam, E Zahavi и др ) При этом контактная задача при конечных деформациях по части постановки обычно рассматривалась либо как задача на условный экстремум некого функционала (для упругих материалов), либо как квазивариационное неравенство скоростного типа (в общем случае) Существуют также невариационные постановки, но к ним прибегают сравнительно редко

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

Важно учитывать еще два обстоятельства Во-первых, при работе с фирменными программными продуктами предлагаемые пользователю способы описания движения абсолютно жестких тел могут оказаться непригодными в каком-то частном случае Например, пакеты ANSYS и NASTRAN использует декартовы координаты выделенной точки тела ("pilot node") и углы поворота вокруг нее, что не позволяет решать некоторые задачи со сложным нагружени-ем Поэтому желательно, чтобы пользователь мог самостоятельно задавать системы обобщенных координат для абсолютно жестких тел Во-вторых, в последнее время начал проявляться интерес к микроструктуре деформируемых тел

--- — —- - - --- - " Это предъявляет до-

полнительное требование к контактным алгоршмам, а именно высокую точность расчета в тонком приконтактном слое при конечных деформациях, по-

скольку в технологических процессах зачастую именно там формируется микроструктура материала

Таким образом, разработка более экономичных и простых в программировании способов решения остается актуальной проблемой

Работа выполнялась в соответствии с научно-технической программой "Высокие технологии высшей школы" (утверждена приказом № 486 от 20 03 96 "Об утверждении перечня минвузовских научно-технических программ на 1996 г "), проектами "Исследование пластического течения металла при локальном и комплексном нагружении" и "Исследование характера пластического течения металла при получении тонкостенных осесимметричных деталей методом валковой штамповки", выигравшими конкурсы грантов в 1996 и 2000 г г соответственно, проектом "Исследование пластического течения металла при изготовлении деталей методом валковой штамповки", вошедшим в разовый заказ-наряд в 1999 г, проектом "Исследование напряженно-деформированного состояния и характера пластического течения металла в разделительных и формообразующих операциях при локальном деформировании", вошедшим в единый заказ-наряд в 2000 г , проектом "Исследование пластического формоизменения при комплексном локальном нагружении объекта деформирования", выполненном по заданию Министерства образования на проведение научных исследований в 2005-2007 гг

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

1 Получить скоростной вариационный принцип квазистатического равновесия системы контактирующих тел, явно включающий в себя не только тензорные поля, определяющие напряженно-деформированное состояние (НДС), но также и обобщенные координаты и силы для абсолютно жестких тел

2 Построить модель описания процесса контактного взаимодействия на основе предложенной формы вариационного принципа и системы эволюционных уравнений

3 Выполнить пространственную дискретизацию полученной модели на основе конечноэлементной аппроксимации Разработать методику решения полученной системы обыкновенных дифференциальных уравнений на основе метода Рунге-Кутта

4 Создать и отладить программу для численного решения контактных задач

5 На серии численных тестов продемонстрировать достоверность предлагаемого метода и его преимущества перед существующими

Научная новизна

1 Предложена модификация вариационного принципа скоростного типа, позволяющая описать процесс взаимодействия жестких тел с телами, испытывающими произвольные деформации

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

3 Получено уравнение эволюции границы поверхности контакта Методы исследования:

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

Кроме этого, для сравнения использовался расчет методом множителей Лагранжа с добавками (augmented Lagrangian method), выполненный коммерческим пакетом ANSYS, показавший результаты, сходные с полученными излагаемым методом

Достоверность результатов обеспечивается сравнением с точными решениями двух тестовых задач Герца, а также с решением третьей тестовой задачи, полученным одним из классических методов, поддерживаемым коммерческим пакетом ANSYS

Практическая значимость и реализация работы

Разработан и реализован как пакет программ на языке FORTRAN 95 метод численного решения контактных задач с трением, ориентированный на задачи обработки металлов давлением, в которых одно упругопластическое тело взаимодействует с несколькими абсолютно жесткими телами Метод может оказаться полезным и в других областях, например, в механике соударения твердых тел На защиту выносится

постановка и метод решения задачи стационарного взаимодействия абсолютно жестких тел с конечно-деформируемым упругопластическим телом

Апробация По содержанию диссертационной работы были сделаны доклады на международном научно-техническом симпозиуме "Механика и технология в процессах формоизменения с локальным очагом пластической деформации", октябрь 1997, ОрелГТУ, на международной научно-технической конференции "Ресурсосберегающие технологии, оборудование и автоматизация штамповочного производства", ноябрь 1999, ТулГУ, на международной научно-технической конференции "Информационные технологии в обработке давлением", апрель 2008, Украина, Краматорск Донецкой обл , на ежегодных научно-технических конференциях в ОрелГТУ

Публикации. По теме диссертационной работы опубликованы б печатных работ, среди которых 2 статьи в центральных научных рецензируемых изданиях, 3 статьи в различных межвузовских сборниках научно-технических трудов, 1 тезисы доклада на научно-технической конференции

Структура и объем работы Диссертационная работа изложена на 97 страницах машинописного текста и содержит 15 рисунков Состоит из введения, трех глав, заключения и общих выводов по работе, списка используемых

источников, включающего 59 наименований работ отечественных и зарубежных авторов.

ОСНОВНОЕ СОДЕРЖАНИЕ РАБОТЫ

Во введении обоснована актуальность темы, сформулирована цель работы, научная новизна, методы исследования, определена практическая ценность работы, приводятся данные о публикациях, структуре и объеме диссертации Сюда помещена также и обзорная часть, состоящая из анализа известных методов решения контактных задач Отмечается вклад отечественных ученых А С Кравчука, В С Давыдова, Е Н Чумаченко, Э Р Гольника, Н И Гундоро-вой, А А Успехова Обзор не ставит целью полноту охвата литературы по контактным задачам (такие публикации уже существуют) Вместо этого на примере нескольких статей из современной научной периодики показаны характерные особенности каждого из трех известных численных методов, называемых автором диссертации "классическими" метода штрафа, метода множителей Лагранжа, метода множителей Лагранжа с добавками (реализованы, в частности, в хорошо известных коммерческих пакетах АЫБУБ и NASTRAN) Показаны их преимущества и недостатки, заключающиеся в следующем

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

Достоинствами метода штрафа являются относительная простота как теоретических основ, так и программирования, а также применимость к широкому кругу задач Кроме этого, как указывается в документации к пакету программ АШУБ, при сильном искажении конечноэлементной сетки данный метод может превосходить своего ближайшего конкурента - метод множителей Лагранжа с добавками

Метод множителей Лагранжа Условие непроникновения выполняется точно, однако этот метод приводит к матрицам систем разрешающих уравне-

ний, обладающим одновременно несколькими нежелательными свойствами знаконеопределеиность, несимметричность, плохая обусловленность, присутствие нулей на главной диагонали Это нередко приводит к сравнительно большому объему вычислений и медленной сходимости, а наличие нулей на главной диагонали затрудняет применение итерационных методов решения Отмеченные недостатки делают метод множителей Лагранжа неконкурентоспособным при конечных деформациях

Достоинством метода является отсутствие дополнительной переменной g (зазора), описывающей несуществующую степень свободы и, как следствие, отсутствие подгоночных параметров вроде штрафных множителей и точное выполнение условия непроникновения Поэтому ошибка численного решения в тонком приконтактном слое сравнительно низка При малых деформациях, когда время счета не столь критично, метод множителей Лагранжа превосходит остальные Превосходство еще ощутимее, если к точности решения таких задач предъявляются повышенные требования

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

За критическим обзором известных подходов к решению контактных задач следует краткое описание предлагаемого метода с перечислением его свойств Достижение этих свойств и составляло задачи исследования

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

С этой целью вводятся различные системы координат, необходимые для описания кинематики контактирующих тел

1) неподвижная (вообще говоря, криволинейная), общая для всех тел, обозначим через у = координаты точек деформируемого тела, а через

У5(1) = ^«(О'^СО'^^О)) ~~ координаты точек абсолютно жесткого тела в этой

системе, причем вариационные принципы даются в ковариантном виде, а ко-нечноэлементная дискретизация выполнена в декартовых координатах),

2) подвижная лагранжева, связанная с деформируемым телом (необходима для определения тензора деформации),

3) система обобщенных координат ,<76(;) для п абсолютно жестких

тел (1 < 1 < п), определяющих положение этих тел в пространстве,

1 2

4) подвижные криволинейные £ (;),£ (,) на поверхностях и абсолютно жестких тел (1 < 1 < п), жестко связанные с этими телами, так что координаты материальных точек поверхностей жестких тел в этих системах остаются неиз-

1 ? з

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

£ > 0 - условие непроникновения,

£ = 0 - условие контактирования,

3 12 £ =0,£ = 0-условия прилипания,

3 1 2

£ = 0, £ ^ 0 или £ * 0 - условия скольжения

Кроме того, для точек поверхности каждого абсолютно жесткого тела верно

Общее число систем координат равно 2+2«, где п - число абсолютно жестких тел Таким образом, даже в простейшем случае п=1 необходимы 4 системы координат

Затем даются необходимые понятия кинематики деформируемых сред Пусть г - радиус-вектор произвольной точки неподвижной глобальной системы

Эг I 1 2 3 \ отсчета, э, =-(1 < г < 3) - локальный базис системы координат Iу ,у ,у I,

ду'

э' (1 < 1 < 3) - двойственный ему базис в пространстве ковекторов, так что Эу э' = Тогда градиентом векторного (или ковекторного) поля а(у) называется тензор второго ранга Уа, равный

Да

Уа = э'-= У,ауэ! ®э, = У.а.э' 0Э7 (1)

а ' У ' / 4 '

ду

В чскаркжоп оркн опальной спсчеме коордиил (y',v~,a') е базисом

(l4 Кт,К-;) фаДПСШ IIMCCI вид

Dx'

причем iici различия между верхними и нижними индексами Введем 1акже ipaiiLiioiiiipoüaiinuii ipannein Va' , равнин

V.

i1 = )' = V,a/3/ О / =V,Í,/)/ ® >' (2)

( г

В декарюьои opioi опальной lucíеме коордшш он выражаекя как

Va7 =—!~k'Qk', дх1

причем снова nei различия между верхними и нижними индексами Заметим,

чю ich,ору Va' i, оючествеипой литературе cooTßeiußyei leinop Va в зарубежной и паоборо!

Дсформацпеп скороеш называйся /епзор вюрого ранга

<l=~(Vv + VvT) , (3)

1де ^ ееп> поле скорости макрпальных точек сформируемого гела Вихрем на-зываекя кпзор вюрого ранга

w=|(Vv' -Vv) (4)

Да нес даекя вывоа нового скоростного вариационного принципа квази-eiai пчееко! о равновесия системы коптакшрующич 1ел Сначала контактирующие 1ела (деформируемое и абсолютно /кесткое) рассматриваются по отдельное in, и чля кажцого записывается известный из литературы вариацион-нып принцип равновесия Для деформируемо! о тела ло - скоростной принцип, впервые полученный в 1983 г О Л Толокоппиковым, А А Маркиным н ВФ Астаповым и впоак истин видоизмененный А А Поздеевым, П В Трусовым и 10 II Пяшпиым (изменение заключается в выражении материальной производной через коро1ацпопп)ю)

jj|Y>Vs (' <v — el a-a w + w с~а w + (V v)a)dV =

'«" , (5)

= ^ p(V v-n cl n))dS

роцс i Bei ni un принципу Журдепа Пепзвесшым здесь являемся поле скороыи

ve wV ( wV

функциональное прост рапс изо Соболева), iiinei рпрование

проп ¡вочи 1ся по ак |уал1,пым обьему тела У{1) и площади кошакта S(t), соответ-ииующпм момеш) времени / a,tl,w - leibopbi неппшого напряжения Кошп, де(|)ормацп11 скорое i п п вихря р, п — давлешн и внешняя нормаль, о i носящиеся к ак lyajii поп моверхпос ni KOinaKia S(l), i - короицноппая производная

в определяющем соотношении для материала деформируемого тела, связанная со спином £1 = (9 + w (в частности, если используется производная Яумана, то w=0) В качестве примера физически корректного подхода к построению определяющих соотношений приводятся работы проф ТулГУ А А Маркина Условия варьирования 1) öv = 0 на поверхности закрепления Su (t), 2) на поверхностях прилипания Sw(t) и скольжения Sa(t) выполнены условия прилипания и скольжения, которые в явном и простом виде приводятся ниже

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

Равновесие абсолютно жесткого тела описывается принципом виртуальной работы, которому автор придает скоростную форму

Q Sq+ {{¿v, (ps+pä(V v-n d n))dS = 0 , (6)

m

где q = |у1,. ,<76) - обобщенные координаты абсолютно жесткого тела (считаем для краткости, что такое тело одно, общий случай очевиден и сводится к расширению размера тензоров q и Q до 6ri), причем для некоторых из них задан закон перемещения (пусть L — множество их номеров), для остальных - закон нагружения (/={1, ,6}\L - множество их номеров), Q = (ßi, - энер-

гетически сопряженные им обобщенные силы, p s - внешнее по отношению к абсолютно жесткому телу давление Ясно, что обобщенные скорости с номерами к е L не варьируются, поэтому соответствующие слагаемые в первой свертке (6) обращаются в ноль

С помощью условий сопряжения в зоне контакта (силового р = -ps, выражающего третий закон Ньютона, а также кинематического v = \t + \s, где vt есть скорость скольжения) вариационные принципы (5) и (6) можно объединить в один новый принцип

j[J<5Vv (ro-d <т-<т <a + w о-д w + (V \)a)dV = Q c5q + V{t)

+ \\övt (p + p(V v-n d nj)dS (7)

Sit)

В координатной записи последний интеграл, как это показано в тексте диссертации, имеет вид

II &

5(0

а , +

" / ' . 7 < гЛ\ . "

j.*,(ajmg + bjr£) + ajg —

dS

где = g„| =е„ eJ - метрика на по-

верхности абсолютно жесткого тела (ввиду недеформируемости тела она не за-

1 2

висит от /), с 1,02 - касательный базис системы координат Е, ,

12 3

ненты ковектора давления в базисе е ,е ,е , тензоры а и Ь определяются соотношениями с-1 -/¡е' и у!

Условия варьирования переменных следующие

1) 5ук = 0 на поверхности Би{(), где \ <к<Ъ,

2) = 0 на поверхности где 1 < г < 2,

3) на поверхности Ба (?) вариация скорости скольжения <5£'е, параллельна касательной составляющей давления Е,с', где 1 < г < 2,

4) &11 =0, где 1еЬ

Заданными (известными) переменными являются

1) скорости обобщенных сил (')> где / е /,

2) обобщенные скорости где / е Ь,

3) скорости .у*(г) точек поверхности 5и(/), где 1<&<3

Напомним, что [у1 ,у2, у3 ] - глобальные криволинейные координаты, которые, в частности, могут быть декартовыми Распределенные неконтактные нагрузки пока не рассматриваются, чтобы не загромождать формулы Они будут учтены ниже по стандартной методике МКЭ

Помимо вариационного принципа (7), выражающего условие квазистационарного равновесия, в постановку задачи необходимо включить кинематиче-

т

ские выражения тензоров Уу, Уу , <1, со, \у через поле скорости, а также определяющее соотношение, уравнение эволюции тензора напряжения, начальные и граничные условия Все необходимые кинематические выражения выписаны в начале раздела за исключением тензора \у В приводимых ниже численных примерах полагалось \у=0, что означает использование производной Яумана Определяющее соотношение запишем в виде

= Ъ (8)

Оно является достаточно общим и позволяет учесть все подходы к построению таких соотношений В частности, уравнение га = Б (1 легко приводится к указанному виду, если учесть, что с! = +

Дадим конкретизацию (8) для приводимых ниже численных примеров В них используется соотношение ' <т = Б 41 с производной Яумана Тензор Б в глобальных декартовых координатах имеет вид

3

ЕЧ*1

ст2(1 + Я'/(ЗС))

при активном нагружении

Е1,к1 при обратимом нагружении

где Е1^ — изотропный тензор упругости, девиатор тензора напряжения, а — интенсивность напряжения Условие перехода из упругого состояния в пластическое а=Н(£р), где II{¿р) - кривая упрочнения, ¿;р - параметр Удкви-

ста Поскольку активное нагружение сопровождается упругим деформированием, приведем правило разложение деформации на упругую и пластическую составляющие Именно упругая деформация считается малой и принимается аддитивное разложение скорости

с1 =йе+Ар (9)

на упругую Ае и пластическую йр Для упругой деформации при активном на-гружении принимается закон

^ = Е йе (10)

Условие активного нагружения о с! > 0, (11)

условие начала разгрузки а с1р<0 (12)

Таким образом, для численных примеров берется теория изотропного упрочняющегося материала с малой упругой деформацией

Эволюционное уравнение напряженного состояния имеет вид

а=;а-а £1 + П а (13)

Начальные условия, помимо исходного положения тел, формулируем как нена-груженность в момент времени (=0

°|/=0= 0, <4,0= 0 (14)

Наконец, граничные условия в контактных задачах задаются неявно через закон трения В численных примерах принимался изотропный закон Кулона с постоянным коэффициентом к

|Р,МР„| (15)

Итак, постановка вариационной задачи полностью завершена Полученная система уравнений (1)-(4), (7)-(15) описывает процесс квазистатического взаимодействия жестких тел с конечно-деформируемым упругопластическим телом Отметим две ее особенности Во-первых, одновременно рассматриваются как поля (скорость, напряжение, давление), так и дискретные переменные (обобщенные координаты и обобщенные силы) Во-вторых, обобщенные координаты абсолютно жестких тел и обобщенные силы явно входят в постановку задачи, что существенно упрощает программную реализацию

Аналитическое решение системы (1)-(4), (7)-(15) невозможно, поэтому необходимы численные методы В качестве такового был выбран метод конечных элементов (МКЭ), хорошо зарекомендовавший себя в области нелинейных задач механики твердого тела. Применению МКЭ к решению данной системы уравнений полностью посвящена вторая глава

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

независимыми эффектами "наплывом-откатом" материала деформируемого тела, а также скольжением по поверхности абсолютно жесткого тела

Вторая глава описывает процедуру конечноэлементной пространственной дискретизации системы уравнений (1)-(4), (7)-(15) и методику решения получающейся в результате разрешающей системы уравнений Сюда включена также и программная реализация

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

vl =хг =х'пН"{х),

где 1<г<3,l<n<N, Нп{х) - базисная функция, связанная с и-ым узлом

A„(xj,,x^,x^J КЭ сетки, N- число этих узлов На поверхности контакта S(t) перейдем к определенным ранее криволинейным координатам

где 1 < j < 2, ) _ криволинейные координаты узла Ап е S(t) Правая часть

(7) после этого принимает вид

kel S(t)

+ Sf¿qr \\aljrE¡HndS

3 1

dS +

(16)

so)

где l</<3,l<r<6,l<i,j<2,l<íw,M<jV,

S] = \\H»~,JS S(l)

есть скорости узловых сил трения Следующим шагом дискретизируем левую часть (7), используя определяющее уравнение (8) Имеем

Ш&'7 iD'ükx\ - {djk - wjk)ah - о* (щ, + wh))dV =

т

= = (17)

т

где \<i,j,k,l <3,1<m,n<N, запятой обозначено дифференцирование по декартовым координатам, В^ и С"^ - коэффициенты пропорциональности, возникающие при дискретизации тензоров d - w и ш + w

dlk~wjk= xnB"jk > °>ki + wh = xnC¡Ja '

в частности, если в определяющем уравнении используется коротационная производная Яумана, то \у=0, 2В"/к = + %//" , 2С"ь = - 8иИ\ ,

где символ Кронекера,

Кшп = до нт{в,,1кнп _ вп 1г - ^С^Г (18)

У( О

есть глобальная матрица жесткости МКЭ Далее показано, как последнюю сумму в (17) представить в форме, включающей обобщенные координаты абсолютно жесткого тела Переходя к одноиндексной нумерации всех степеней свободы д', для узлов Лп 5, Для узлов Ап е 5, выводим из дискретной формы вариационного принципа (7) типичную для МКЭ линейную систему ОДУ

Кт„и"=Рт 09)

где и" - степени свободы в одноиндексной нумерации, ^ - энергетически сопряженные им силы После некоторых преобразований этой системы (заполнения нулями строк матрицы К, соответствующих неварьируемым степеням свободы за исключением элементов этих строк, находящихся на главной диагонали, которые замещаются единицами, а так же заменой в правой части скоростей сил на скорости координат для неварьируемых степеней свободы) получаем разрешающую систему ОДУ Заметим, что силы трения задаются неявно законом трения, однако при вычислениях учесть их не составляет труда Распределенные неконтактные нагрузки (силу тяжести, гидростатическое давление и др) можно внести непосредственно в правую часть (19) согласно формулам, приводимым в учебниках по МКЭ

Приведенная система ОДУ интегрируется на основе многомерного варианта метода Рунге-Кутта 4-го порядка

^ = хп+\ ~хп

Уп+\ = Уп + \ (/о + 2/1 + 2/2 + /з )И /О = /(* п'Уп)

Ь=Ахп+Н>Уп+/2Н)

Выбор именно этой схемы объясняется следующими соображениями Особенностью задач обработки металлов давлением является то, что картина деформирования гораздо сильнее зависит от параметра нагружения, чем от пространственных координат Так, если для КЭ-аппроксимации по пространственным координатам обычно хватает -100 конечных элементов, приходящихся на характерный размер, то число шагов нагружения обычно на несколько порядков

выше, а при линейной аппроксимации по времени на шаге вообще может оказаться неприемлемо большим Отсюда ясно, что процедура интегрирования по времени должна иметь больший порядок точности, чем по пространственным координатам Отметим также, что при решении контактных задач ввиду очень сложных граничных условий и различных нелинейностей могут возникать проблемы с устойчивостью вычислений Следовательно, процедура интегрирования по времени системы разрешающих уравнений должна быть очень устойчивой и допускать широкое варьирование шага интегрирования Этим условиям удовлетворяет метод Рунге-Кутта Наиболее эффективной по соотношению "точность-объем вычислений" считается вышеприведенная схема 4-го порядка При ее использовании в целях упрощения вычислений используется рекомендуемое рядом авторов разделение нелинейностей задачи на физические и геометрические с целью отбросить последние, что может быть оправдано, если шаги нагружения в пошаговом анализе выбраны достаточно малыми, и отсутствует эффект изгибания

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

Третья глава посвящена численным тестам, демонстрирующим достовер-

рис 3 Задача Герца рис 4 Решение задачи Герца

Пример 1 Контактная задача Герца

Классическим тестом для контактных пакетов программ является задача Герца, поскольку она имеет аналитическое решение Ее особенностью являются малые деформации и малые повороты, следовательно, можно использовать

практически любые определяющие соотношения Бесконечно длинная толстая труба внешнего радиуса 0,16 м и внутреннего 0,1 м сжимается двумя параллельными абсолютно жесткими плитами с силои F (см рис 3) Трение отсутст-

Q 7

вует Свойства материала упругий изотропный, модуль ЮнгаЕ = 9,8 10 п/м , коэффициент Пуассона ^=0,3 Требуется найти распределение давления по поверхности контакта для двух значений погонной сжимающей силы

2 2 F = 9,8 10 н/м и F = 19,6 10 п/м Задача имеет горизонтальную и вертикальную оси симметрии Решения в предположении плоского деформирования приведены на рис 4, где они сравниваются с аналитическими Использовалась конечноэлементная сетка с переменным шагом из линейных треугольных элементов Давление рассчитывалось в узлах сетки Погрешность решения для данной дискретизации можно признать удовлетворительной Она примерно соответствует погрешности, полученной другими исследователями, использовавшими метод множителей Лагранжа, который в данном случае, ввиду малых деформаций, предпочтительнее из всех классических

_ j

Пример 2 Заготовка в виде сектора кольца с внешним радиусом 3 10 м,

_2

внутренним радиусом 2,5 10 м и углом л-/8 закреплена по внутренней дуге

АВ и обкатывается по внешней дуге CD роликом радиуса 1,5 10~2л* с центром О (см рис 5, на котором показано исходное положение) При этом центр ролика равномерно поворачивается по часовой стрелке на угол яг/16 вокруг центра сектора, к ролику прикладывается линейно возрастающая от нуля сила прижима, направленная к центру сектора, а также линейно возрастающий от нуля тормозной момент относительно центра О Деформация считается плоской Максимальное погонное значение прижимающей силы Fq =7,5 105 н/м,

максимальное погонное значение тормозного момента Mq =2 10 к м/м Ocio 2

тальныеусловия задачи следующие модуль сдвига (7=2,7 10 п/м .коэффициент Пуассона /л =0,32, трение подчиняется закону Кулона с коэффициентом 0,2, пластическое поведение материала считается изотропным с линейным упрочнением, начальный предел текучести ctq = 1,93 10 «/д! , модуль упрочне-8 2

ния//'=9,62 10 н/м Похожие задачи возникают при исследовании новой

технологии поверхностного упрочнения деталей Расчеты производились на последовательно измельчаемых сетках до достижения следующего критерия при измельчении сетки в 2 раза интенсивность напряжения меняется в пределах 10% в той области, где она не слишком мала (не менее 10% от максимальной) Использовались линейные треугольные конечные элементы Результаты расчета на самой мелкой сетке представлены на рис 6а, 7а в Та же задача решалась пакетом ANSYS на сетке

Рис 6а Изолинии интенсивности

О 'У

напряжения [10 н/м ] ниже стд, новый метод 1 -0,43,2-0,87, 3-1,3, 4-1,73

Рис 66 Изолинии интенсивности

о 9

напряжения[10 к/м ] ниже стд, пакет ANS YS 1 -0,4,2-0,8, 3-1,2,4-1,55

Рис 76 Изолинии интенсивности

о 9

напряжения [10 н/м ] выше ctq, пакет ANS YS 1-2,5,2-2,9, 3-3,3,4-3,9

примерно с той же густотой узлов, но с использованием треугольных 6-узловых квадратичных конечных элементов (сетки из линейных треугольных элементов не рекомендованы разработчиками ANSYS - см документацию для ANSYS) Метод решения - множителей Лагранжа с добавками как наиболее эффективный из классических Результаты расчета представлены на рис 66, 76 Перейдем к сравнению результатов Прежде всего отметим, что параметр Удквиста не превысил 0 25, что говорит о приемлемости производной Яумана в данном расчете Для удобства изолинии интенсивности напряжения ниже и выше начального предела текучести показаны на разных картинках и в разных масштабах

Дуга контакта ЕР всюду выделена жирной линией Видно, что изолинии рис 6а, 7а и рис 66, 76, соответствующие близким значениям интенсивности напряжения (эти значения различаются на 9-12%), имеют сходную форму Учитывая вычислительную сложность примера 2, согласование изолиний можно признать удовлетворительным и достаточным для подтверждения работоспособности предлагаемого метода

Пример 3 Цель данного примера - продемонстрировать возможности, недоступные пакету Кольцо переменной толщины внешнего ра-

—2 —? диусаЗ 10 ми внутреннего 2,5 10 м с центрами С \ и С2 внешней и внутренней окружностей, отстоящими на 10 м, обкатывается 4 роликами по симметричной относительно прямой С\ С 2 схеме - см рис 8 Обратим внимание, что С] и С2 - начальные положения центров внешней и внутренней границ кольца, поскольку в дальнейшем кольцо деформируется Центры двух роликов 0\ и С>2 поворачивается вокруг С\ навстречу друг другу, причем 0\ остается на неизменном расстоянии от С) К центру О2 прикладывается

примере 3 [10Й н/м2] 1 - 1,93,2- 1,72,3 -0,87

направленная к С\ прижимающая погонная сила Кроме того, к роликам с центрами 0\ и О2 прикладываются одинаковые по величине и противоположные по направлению тормозные погонные моменты относительно их центров Оставшиеся два ролика движутся и нагружаются симметрично первым двум Точкам заготовки на оси симметрии С\С2 запрещается перемещаться в вертикальном направлении и разрешается в горизонтальном Все остальные условия взяты из примера 2

Заметим, что из-за переменной толщины кольца задачу нельзя переформулировать таким образом, чтобы рассматривать движение центров роликов 0\ и

С>2 в декартовых координатах Следовательно, пакеты ANSYS и NASTRAN, по крайней мере в своем стандартном исполнении, решить пример 3 не в состоянии

Результаты расчета представлены на рис 9 Параметр Удквиста не превысил 0 03, что говорит о приемлемости производной Яумана в данном расчете В силу симметрии показана половина заготовки Зона пластичности ограничена изолинией №1 Видно, что один ролик пластическую деформацию не вызывает, тогда как под воздействием другого она развивается не только под прокатанной поверхностью, но и с обратной стороны заготовки С практической точки зрения это свидетельствовало бы о необходимости изменения данной схемы упрочнения

ОСНОВНЫЕ РЕЗУЛЬТАТЫ И ВЫВОДЫ

1 Решена новая актуальная научная задача, состоящая в разработке и численном исследовании модели контактного взаимодействия абсолютно твердых тел с упругопластическим телом, испытывающем конечные деформации

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

3 На основе конечноэлементной аппроксимации выполнена пространственная дискретизация системы уравнений, описывающей процесс контактного взаимодействия

4 Разработана методика решения полученной системы обыкновенных дифференциальных уравнений на основе метода Рунге-Кутта

5 Получено уравнение эволюции границы поверхности контакта, явно определяющее эту поверхность

6 Создана и отлажена FORTRAN программа для численного решения контактных задач

7 На серии численных тестов продемонстрирована достоверность предлагаемого метода и его преимущества перед существующими К тому же предлагаемый метод выгодно отличается от существующих (штрафных функций (Р), множителей Лагранжа (L), множителей Лагранжа с добавками - augmented La-grangian method (А), используемых в пакетах ANSYS и NASTRAN) эффективностью как при больших, так и при малых деформациях (в отличие от Р, L, А), а также высокой точностью решения в тонком приконтактном слое благодаря выполнению условия непроникновения (в отличие от Р и А)

8 Благодаря эффективности и простоте программирования новый метод может быть полезен исследователям, располагающим даже небольшим бюджетом

ОСНОВНЫЕ РЕЗУЛЬТАТЫ ДИССЕРТАЦИИ ОПУБЛИКОВАНЫ В РАБОТАХ

1 Морев П Г Конечноэлементный анализ упругопластических моделей в процессах ОМД//Сборник научных трудов ОрелГТУ 1996, Т 9, С 42-47

2 Морев П Г Использование схемы Рунге-Кутта в конечноэлементном анализе процессов ОМД// Исследования в области теории, технологии и оборудования обработки металлов давлением Межвузовский сб научных трудов Орел ОрелГТУ-ТулГУ,1998, С 134-139

3 Морев П. Г. Учет контакта инструмент-заготовка в конечноэлемент-ных моделях// Известия ТулГУ, сер. Механика деформируемого твердого тела и обработка металлов давлением. 2003, вып. 1, С. 51-56.

4 Море» П. Г. Вариант метода конечных элементов для контактных задач с трением// Известия РАН, сер. Механика твердого тела. 2007, №4, С. 168-182.

5 Голенков В.А, Морев П Г , Радченко С Ю Методы математического моделирования и новые задачи ОМД// Тематический сборник научных трудов "Совершенствование процессов и оборудования обработки давлением в металлургии и машиностроении" 2008, Украина, г Краматорск Донецкой обл , ДГМА, С 15-19

Подписано к печати 27 05 2008 г Формат 60x84 1/16 Печать офсетная Объем 1,0 уел пл Тираж 100 экз Заказ № 1184

Отпечатано с готового оригипал-макета па полиграфическом базе Орловского государственного технического университета 302020, г Орел, Лау горе кое шоссе, 29

Введение.

1. Постановка вариационной задачи.

1.1 Описание рассматриваемого класса задач.

1.2 Постановка контактной задачи в общем виде.

1.3 Кинематика и определяющие соотношения конечного упругопласти-ческого деформирования.

1.4 Вариационный принцип и полная система уравнений.

1.5 Учёт изменения локального базиса.

1.6 Зависимость поверхности контакта от времени.

27 Пространственно-временная дискретизация.52,

2.1 Пространственная конечноэлементная дискретизация.

2.2 Вычисление узловых и интегральных (обобщённых) сил.

2.3 Временная дискретизация на основе метода Рунге-Кутта и интегрирование системы разрешающих уравнений.

2.4 Программирование и общая стратегия расчёта.

2.5 Описание алгоритма.*.j.:.

3. Численные примеры.

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

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

2)Сложная, меняющаяся по времени картина участвующих в расчёте тензорных полей заставляет прибегать к сильно неравномерной дискретизации, сгущающейся в областях большого градиента этих полей. Ввиду остаточной пласти- -: ческой деформации, меняющей свойства материала, дискретизацию-приходится измельчать не только в области контакта тел в текущий момент времени, но и вблизи участков поверхности, находившихся в зоне контакта до этого момента. В результате дискретизация может оказаться настолько мелкой, что достижение погрешности в 10% на персональном компьютере за разумное время будет пробле- -матичным уже для 2-мерных задач. ~ , - •

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

Ввиду своей важности и сложности контактные задачи привлекали большое число исследователей как в нашей стране (А.С. Кравчук, B.C. Давыдов, Е.Н. Чу-маченко, Э.Р. Гольник, Н.И. Гундорова, А.А. Успехов и др.), так и за рубежом (G.

Pietrzak, A. Curnier, F. Armero, E. Petoch, P. Alart, M. Barboteu, F. Lebon, D. Bar-lam, E. Zahavi и др.). В результате было разработано немало алгоритмов и пакетов программ (например, [1-7], пакеты ANSYS, NASTRAN, DEFORM). В принципе, они охватывают почти весь диапазон технологических задач, однако на практике получение достоверных результатов для сложных процессов деформирования может оказаться серьёзной проблемой, связанной с подгонкой нескольких параметров, отсутствующих в постановке задачи и влияющих только на вычислительный процесс. Эти параметры могут меняться от задачи к задаче, поэтому их приходится каждый раз подгонять заново. При этом обычным делом является расходимость пошагового расчёта после довольно большого числа шагов, никак не связанная с потерей устойчивости реальной механической системы контактирующих тел. Кроме этого, предлагаемые алгоритмы чрезвычайно сложны, и для своей программной реализации требуют участия как квалифицированных математиков, так и квалифицированных программистов, что может оказаться неприемлемым для исследователей, располагающих небольшим бюджетом.

Таким образом, разработка более экономичных и простых в программировании способов решения остаётся актуальной проблемой. В диссертации предложен новый подход для важного частного случая формоизменения упругопластиче-ского тела при контактировании с абсолютно жёсткими телами, превосходящий известные методы решения сразу по нескольким параметрам. Подход использует метод конечных элемента (МКЭ) и позволяет учитывать сложные законы движения и нагружения абсолютно жёстких тел, а так же большие упругопластические деформации. Причём нагрузки могут быть не только распределёнными, но и интегральными (т. е. суммарными) для каждого из абсолютно жёстких тел. Прежде чем перейти к его изложению, остановимся сначала на наиболее популярных существующих методах. Их три: метод штрафных функций, метод множителей Jla-гранжа и метод множителей Лагранжа с добавками (augmented Lagrangian method). Нет необходимости делать представительный обзор по практической реализации этих методов, поскольку такие обзоры уже существуют - см. [52] и приводимые там ссылки. Вместо этого вкратце изложимшекоторыё^современные работы,; характеризующие все три указанных подхода.

Метод штрафных функций. Его особенностьюявляется использование до- . полнительной переменной, характеризующей пару контактирующих тел, - зазора между ними; Дадим определение этой переменной, следуя [5]. В данной работе использован лагранжев подход. Пусть X и х(/) = — начальное и текущее положения произвольной материальной точки одного из двух контактирующих тел, Z h z(t) = v.(2) (Z,t) - то же самое для другого, Г^ и Г^ - границы этих тел, привязанные к отсчётной конфигурации: Определим проектирующее отображение точек границы одного тела на границу другого Y(Y): Г^ -» Г^ как

Y(X,0 = arg min ш(1)(Х,0 - V(2)(V,0 , , - . . УбГ<2> \ , . , где X e Г^, а верхний индекс указывает на одно из двух контактирующих тел. •: -Зазор^междутелами! V(1)(X,0-V(2)(V(X,0,0 = ,§rn , • - - - ' ' где n— внешняя нормаль к поверхности контакта (г^"^,ifj в.текущей конфи- -гу рации. Необходимо отметить, что данный метод не предполагает.выполнение г,, условия непроникновения, откуда, собственно; и появляется дополнительная пёг ременная g. Вследствие этого поверхности контакта текущей конфигурации не обязаны совпадать. Корректность приведённого определения зазора g следует из того; что, по определению отображения X -> Y(X, t) , вектор в левой части последнего равенства действительно направлен по нормали к поверхности в точке Y(X,/).

Далее авторами [5] формулируется вариационный принцип для системы контактирующих тел, в котором варьируется, в том числе, и переменная g. Соответствующая этой переменной виртуальная работа является вариацией некого функционала, называемого штрафным, входящего в функционал задачи отдельным слагаемым и принимающего нулевое значение, если условие непроникновения в точности выполняется и очень большое значение, если это условие хотя бы слегка нарушается. Тогда решение вариационной контактной задачи, обеспечивающее экстремум функционала задачи, не может сильно отклоняться от условия непроникновения. Необходимая точность расчёта достигается подбором штрафного функционала исходя из конкретных условий. При дискретизации он становится (штрафной) функцией, с чем и связано название метода.

Из сформулированного авторами [5] вариационного принципа выводится схема пошагового численного интегрирования. Для нахождения зазора g используется дополнительная схема интегрирования, основанная на том, что, как можно показать, g = (v^(x,o-v(2)(Y(x,o,o)-n , где X е Г^, V^ - поле скорости поверхности /-го тела, привязанное к отсчёт-ной конфигурации. Пусть g/, g/+i - величина зазора на последовательных шагах нагружения. Тогда нормальная компонента давления рп в зоне контакта, энергетически сопряжённая переменной g, находится по схеме

U(gi+1)-U(gi)

Рп =1

Si+1 - Si r\ if gi+l*gi

U' ~(gi+\+gi) if gi+\=gi > / для неотрицательного регулярного штрафного потенциала U(g)~. В приведённых авторами численных примерах полагается

U(g) = hg2 if g<0

0 if g > о при большой, выбираемой пользователем, величине штрафного параметра к. Эти примеры демонстрируют поведение различных пар соударяющихся упругих тел, сопровождаемое конечным деформированием.

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

Метод мноэ/сителей Лагранжа. Автору известны 3 относительно недавние работы [6, 53, 54], причём две из них ограничиваются частным случаем малых деформаций. Учитывая, что, по мнению большинства специалистов, этот метод уступает по вычислительной эффективности своим конкурентам в случае конечных деформаций, сделаем лишь общие замечания. Если [54] основана на классической версии метода, то [6, 53] привлекает методы математического программирования. А именно: сначала формулируется контактная задача как вариационное неравенство со специфическими условиями варьирования входящих в него полей; затем производится дискретизация и переформулировка в задачу теории оптимизации. При этом условия на поверхности контакта преобразуются в необходимые и достаточные условия Куна-Таккера. Наконец, новая задача решается с помощью известных методов математического программирования, и в результате находится решение исходной задачи.

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

Метод мноэ/сителей Лагранжа с добавками. Этот сравнительно новый подход оказывается существенно сложнее двух предыдущих из-за методов исследования, заимствованных в негладком выпуклом функциональном анализе. Достаточно отчётливое представление о нём даёт работа [55], которой мы и будем следовать. Начнём с того, что изложим основные идеи.

Многие задачи в физике, механике и экономике могут быть сформулированы как минимизация суммы двух функционалов infF(v)=mf(/(v) + ^(v)) , veH v&H где if есть гильбертово пространство-соскалярным-произведением, обозначае- - -. мым далее как (у); /,¥: Н -> R есть выпуклые нижние полунепрерывные функционалы и/непрерывно дифференцируем; R = R и {- оо,+оо} есть расширенная числовая прямая; функция g такова, что сохраняется выпуклость функционала T(v) = 4>(g(y)). В приведённой сумме 1Р обычно выражает ограничения, наложенные на v или на/ Рассмотрим задачу минимизации с ограничением inf /(v), veC где~С есть непустое замкнутое выпуклое подмножество Н, которую переформулируем в задачу без ограничений inf F(v), где v&H

ГО ifveC F(v) = /(v) + /c(v), /c(v)= .

00 if V £ С

Удобный путь решения новой задачи минимизации - это преобразовать её в задачу минимакса с седловой точкой, которая состоит в нахождении такой пары (и,А), что где Г- другое гильбертово пространство, Y* - сопряжённое ему, X — дуальная переменная, называемая лагранжевым множителем, l{u,A) = -swp{{A,p)-F{u,p)) = inf (F{u,p)-{A,p)) -peY P^Y есть лагранжиан задачи, . . F(u,pj=f(uj+4>(g(uj+p) есть возмущение функционала F.

Трудности с решением задачи поиска минимакса лагранжиана возникают, как только оказывается недифференцируемым. Один из самых популярных и эффективных методов обойти эту трудность основан на так называемой технике лагранжиана с добавкой. Ключевая идея состоит в замене 1(и,Л) на регуляризованный функционал L(u, А) с седловой точкой, называемый лагранжианом с добавкой, который сохраняет все существенные свойства /(и, Л) и может исследоваться стандартными методами. Рассматриваются только два практически важных случая: 1) 4{g{u)) = Ic{g{u)), 2) 4!(g(u)) = Ic(g(u));'v%e Ic есть функционал, сопряжённый Iq (С - выпуклое замкнутое множество). Заметим, что 1) и 2) охватывают все равенства, неравенства и условия типа выпуклости множества, накладываемые как на исходную и, так и на дуальную X переменные. Первый случай представляет минимизацию функционала/с ограничением, наложенным на исходную переменную иеС, где С = {и | g(u)> 0}: inf f(u) = mf(f(u) + Ic(u)) . и&С и

Соответствующий лагранжиан с добавкой имеет следующую форму: lp lp л

9 , \ 2 где disfg{z)= inf || z - q || есть квадрат расстояния между точкой z и выпуклым qeB множеством В, ар есть подстраиваемый параметр регуляризации. Второй случай, соответствующий минимизации с ограничением Я е С, наложенным на дуальную переменную, , inf f(u) = inf (J (и) + 1*с 00) ,

ЛеС и приводит к следующей форме лагранжиана с добавкой:

L(u^) = f(u) + {*,g(u))+£\\gm2~distl& + pg(u)) . . .

2 2 р

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

Таковы основные идеи этого метода. Применительно к контактным задачам он конкретизируется следующим образом. Снова используем лагранжев подход и запишем условие непроникновения с помощью введённых ранее обозначений:

• • g(X,0>0, pn{X,t)>Q, g(X,t)pn(X,t) = 0 , где ХеГ®.' Как оказывается, это условие можно переформулировать на языке функционального анализа следующим образом: нормальная компонента давленияpn(X,t) для X g Г^ является субградиентом недифференцируемого выпуклого потенциала Ir (g(X,t)). Или, что равносильно: зазор g(X,t) для X е Г^ является субградиентом недифференцируемого выпуклого потенциала (pn(X,t)).

Другими ограничениями, требующими учёта, являются условие скольжения и закон трения.Условие скольжения:

Ivt Сх, 0||рг (х, 0 = ||рг (х, t)\\t (х, t), где vt(X,t) и р^(X,/) - скорость скольжения и касательная составляющая давления в точке контакта X е Г^. Закон трения (изотропный): рдх,0|-^(х,0)<0.

В случае закона Кулона получаем k{pn(X,t)) = -fjpn{X,t) с коэффициентом трения ju. Эти ограничения также можно выразить на языке функционального анализа: рДХ,?) является субградиентом недифференцируемого выпуклого квазипо тенциала диссипации ^(уt(X,t)), где С(рп) — двумерное множество касательных компонент давления, допускаемых законом трения, при фиксированном рп. В случае изотропного закона трения С{рп) является кругом. , .

Далее исходная контактная задача формулируется как поиск условного экстремума некого функционала Usys = + и(2) , где слагаемые относятся к контактирующим телам по отдельности. Затем происходит переформулировка в задачу на безусловный экстремум сингулярного функционала

Usys +U& + Я IR Сg(X,t))dT<!> + JJ IciPn)(yt(X,t))drW . r(l) r(l)

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

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

Изложив основные принципы всех трёх методов, перейдём к критическому анализу.

Метод штрафных функций использует подстроечные параметры (один или несколько штрафных множителей), оптимальный выбор которых может вылиться в самостоятельную проблему, требующую экстраполяции результатов нескольких расчётов [8]. Это особенно неудобно для з ад ач, тр ебу ю щи хд лите л ьн ыхко м п ь га-терных вычислений. Ещё одна неприятная : особенность,-- влияние штрафных множителей на обусловленность разрешающей системы уравнений, в результате чего обусловленность может значительно ухудшаться в процессе оптимизации этих множителей. ■ - - - ;; ;, :/

Метод множителей Лагранжа приводит к матрицам, обладающим одновре- . менно несколькими нежелательными свойствами: .знаконеопределённость, несимметричность, плохая обусловленность, присутствие нулей на главной диагонали. Это нередко приводит к'сравнительно большому объёму вычислений и медленной сходимости, а наличие нулей на главной диагонали затрудняет применение итерационных методов решения [7].

Чтобы уменьшить объём вычислений, метод множителей Лагранжа комбинируют с методом штрафных функций и приходят к методу множителей Лагранжа с добавками (augmented Lagrangian method); при этом каждый лагранжев множитель оказывается суммой двух слагаемых, одно из которых пропорционально параметру р регуляризованного лагранжиана. Хотя обусловленность матриц суще- • -ственно улучшается, их всё же рекомендуется подвергать специальному преобразованию с целью дальнейшего улучшения [7]. От метода множителей Лагранжа наследуются также и все остальные нежелательные свойства матриц, а от метода штрафных функций — подстраиваемые параметры.

Необходимо отметить ещё одно нежелательное свойство, присущее всем трём описанным методам. Это - сравнительно невысокая точность в тонком прикон-тактном слое при конечных деформациях. Дело в том, что в последнее время растёт интерес к микромеханике деформируемых сред. Это связано с разработкой новых технологий обработки металлов давлением, когда требуется не только придать детали заданную геометрию, но и управлять её микроструктурой, которая зачастую формируется как раз в приконтактном слое. Если первый и третий методы не обеспечивают высокую точность из-за нарушения условия непроникновения, то второй метод, хотя это условие и соблюдает, оказывается малоэффективным при конечных деформациях, как уже говорилось выше. * Следует упомянуть ещё о двух методах, встречающихся в литературе по кон- . тактным задачам: методе математического программирования и методе вариационных неравенств.

Метод математического программирования имеет отношение скорее к постановке задачи, чем к решению и сводит её к задаче оптимизации с необходимыми и достаточными условиями Куна-Таккера, которая чаще всего решается методом множителей Лагранжа с различными вариациями разрешающего алгоритма [6, 53].

Метод вариационных неравенств, как и метод математического программирования, используется для постановки контактной задачи и не имеет отношения к её решению. Решаются же вариационные неравенства одним из трёх описанных методов.

Метод, предлагаемый в настоящей диссертации, не содержит подгоночные параметры и приводит к хорошо обусловленным матрицам, — их обусловленность обычно не намного хуже, чем у глобальной матрицы жёсткости, получающейся сразу после поэлементной сборки (в противном случае можно выполнить элемен- ■ тарную процедуру масштабирования, т.е. удобного выбора единиц измерения -обобщённых координат). Диагональные элементы оказываются максимальными по абсолютной величине в своих строках или близкими к максимальным. Более того, в упругопластических расчётах на основе формулировок, приводящих к матрицам, представляющим собой сумму положительно определённого симметричного слагаемого и несимметричной добавки, существует практика отбрасывания этой добавки с сохранением сходимости итераций Ньютона-Рафсона, по-видимому, для большинства технических задач. Данный приём позволяет существенно упростить программное обеспечение, используя простые, эффективные и общедоступные алгоритмы [9, 10]. При этом увеличение числа итераций из-за огрубления матриц может компенсироваться увеличением скорости их прохождения за счёт более эффективных матричных алгоритмов (вопрос о количественных оценках на этот счёт требует отдельного исследования и здесь не рассматрива-1 ется). Такая практика совместима-с предлагаемым, методом;, но не-совместима, ; -" например, с методом множителей Лагранжа, включая его.версию с.добавками. \ . ■ . Важно учитывать ещё одно обстоятельство. При работе с фирменными программными продуктами предлагаемые пользователю способы описания движения абсолютно жёстких тел могут оказаться неудобными в каком-то частном слу- • чае. Например, пакет ANSYS использует декартовы координаты выделенной точки тела ("pilot node") и углы поворота вокругнеё. В конце статьи приводится -. • пример, когда эти координаты не подходят. Предлагаемый метод позволяет работать с любыми обобщёнными координатами жёстких тел, причём своими для каждого такого тела и задавать интегральные нагрузки в энергетически сопряжённых обобщённых силах.

Диссертация состоит из введения, трёх глав, заключения и списка литературы.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ И ВЫВОДЫ

1. Решена новая актуальная научная задача, состоящая в разработке и численном исследовании модели контактного взаимодействия абсолютно твёрдых тел с упругопластическим телом, испытывающем конечные деформации.

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

3. На основе конечноэлементной аппроксимации выполнена пространственная дискретизация системы уравнений, описывающей процесс контактного взаимодействия.

4. Разработана методика решения полученной системы обыкновенных дифференциальных уравнений на основе метода Рунге-Кутта.

5. Получено уравнение эволюции границы поверхности контакта, явно определяющее эту поверхность.

6. Создана и отлажена FORTRAN программа для численного решения контактных задач.

7. На серии численных тестов продемонстрирована достоверность предлагаемого метода и его преимущества перед существующими. К тому же предлагаемый метод выгодно отличается от существующих (штрафных функций (Р), множителей Лагранжа (L), множителей Лагранжа с добавками - augmented Lagrangian method (А), используемых в пакетах ANSYS и NASTRAN) эффективностью как при больших, так и при малых деформациях (в отличие от Р, L, А), а также высокой точностью решения в тонком приконтактном слое благодаря выполнению условия непроникновения (в отличие от Р и А).

8. Благодаря эффективности и простоте программирования новый метод может быть полезен исследователям, располагающим даже небольшим бюджетом.

1. Гольник Э.Р., Гундорова Н.И., Успехов А.А. Инкрементальное моделированиефрикционных контактных систем на базе неинкрементального алгоритма// Изв. вузов. Машиностроение. 2000, №3, с. 9-14.

2. Кравчук А.С. Вариационные и квазивариационные неравенства в механике.1. М.: МГАПИ, 1997, 340 с.

3. Давыдов B.C., Чумаченко Е.Н. Метод реализации модели контактного взаимодействия в МКЭ при решении задач о формоизменении сплошных сред// Известия РАН, сер. Механика твёрдого тела. 2000, № 4, с. 53-63.

4. Peter W., Christensen. A. A semi-smooth Newton method for elasto-plastic contactproblems// International Journal of Solids and Structures. 2002. v.39. p. 23232341.

5. P. Alart, M. Barboteu, F. Lebon Solution of frictional contact problems by an EBE ~~ ~" preconditioned/ Computational Mechanics. 1997-. v.20. p. 370-378.

6. D. Barlam, E. Zahavi The reliability of solution in contact problems// Computersand Structures. 1999. v.70. p. 35-45.

7. Джордж А., ЛюДж. Численное решение больших разреженных систем уравнений. М.: Мир, 1984, 333 с.

8. Писсанецки С. Технология разреженных матриц. М.:Мир, 1992.

9. Поздеев А.А., Трусов П.В., Няшин Ю.И. Большие упругопластические деформации: теория, алгоритмы, приложения. М.: Наука, 1986. 231 с.

10. Унксов Е. П., Джонсон У., Колмогоров В. Л. и др. Теория пластических деформаций металлов. М.: Машиностроение, 1983. 598 с.

11. В. Schieck, Н. Stumpf The appropriate corotational rate, exact formula for plasticspin and constitutive model for finite elastoplasticity // International Journal of Solids and Structures, 1995, v. 32, № 24, p. 3643-3667.

12. J. R. Barber, M. Ciavarelly Contact mechanics// International Journal of Solidsand Structures, 2000, v. 37, p. 29-34

13. Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. М: Наука, 1977, 832 с.

14. Кунин С. Вычислительная физика. М. : Мир, 1992, 518 с.

15. Васидзу К. Вариационные методы в теории упругости и пластичности. М. : Мир, 1986

16. Маркин А. А. Об изменении упругих, и пластических свойств .при: конечном -^ ~ деформировании// Известия АН СССР, сер. Механика твёрдого тела, "1990,2, с. 120-126.

17. Левитас В. И. Большие упругопластические деформации материалов при высоких давлениях. Киев: Наукова думка, 1987, 228 с.

18. Svedsen В. On the modelling of anisotropic elastic and inelastic behaviour at large deformation// International Journal of Solids and Structures, 2001, v. 38, p. 95729599.

19. J. Casey, P. M. Naghdi A remark on the use of the decomposition F = Fe • F^ in plasticity// Journal of Applied Mechanics, 1980, v. 47, p. 672-675.

20. Y F. Dafalias The plastic spin// Journal of Applied Mechanics, 1985, v. 52, p. 865871.

21. P. M. Naghdi A critical review of the state of finite plasticity// Zeitschrift fur An-gewandte Mathematik und Physik. 1990, v. 41, p. 315-387.

22. C. Miehe A constitutive frame of elastoplasticity at large strains based on the notation of plastic metric// International Journal of Solids and Structures, 1998, v. 35, p. 3859-3897.

23. C. Miehe A formulation of finite elastoplasticity based on dual со- and contra-variant eigenvector triads normalized with respect to a plastic metric// Computer meth" ods in applied mechanics and engineering: 1998, v.l59, p. 223-260.

24. Jacob Fish, Kamlun Shek Finite deformation plasticity based on the addidive split of the rate of deformation and hiperplasticity// Computer methods in applied mechanics and engineering. 2000, v. 190, p. 75-93.

25. Xiao K, Bruhns O., Meyers A. Logarithmic strain, logarithmic spin and logarithmic J ' rate//Acta Mechanica. 1998~ v. 124, p. 80-105. •

26. Xiao H., Bruhns O., Meyers On objective corotational rates and their defining spin tensors//International Journal of Solids and Structures, 1998, v. 35, p. 4001-4014.

27. Xiao H., Bruhns O., Meyers A self-consistent Eulerian rate type model for finite deformation elastoplasticity with isotropic damage// International Journal of Solids and Structures, 2001, v. 38, p. 657-683.

28. Chung Souk Han, Yangwook Choi, June Lee, R. H. Wagoner A FE formulation for elasto-plastic materials with planar anisotropic yield functions and plastic spin// International Journal of Solids and Structures, 2002, v. 39, p. 5123-5114.

29. A. Menzel, P. Steinmann On the spatiaLformulation of anisotropic multiplicative- -elasto-plasticity// Computer methods in applied mechanics and engineering. 2003, v.192, p. 3431-3470.

30. Ильюшин А. А. Пластичность: Основы общей математической теории. М.: ИПМ АН СССР, 1963, 272 с.

31. Ильюшин А. А. Механика сплошной среды: Учебник для университетов. М.: МГУ, 1978, 287 с.

32. Седов Л. И. Введение в механику сплошной среды. М.: Физматгиз, 1962, 284с.

33. R. Blaheta Displacement decomposition incomplete factorization preconditioning techniques for linear elasticity problems// Numerical and Linear Algebra Application, 1994, v. l,p. 107-128.

34. R. Blaheta, R. Kohut Solution of large nonlinear systems in elasto-plasticity incremental analysys// Journal of Computational and Applied Mathematic, 1995, v. 63, p. 255-264.

35. R. Blaheta, O. Axelson Convergence of inexact Newton-like iteration in incremental -finite element analysis of elasto-plastic problems// Computer methods in applied- • mechanics and engineering, 1997, v. 141, p. 281-295.

36. Черных К. Ф. Введение в анизотропную упругость." M.: Наука, 1998, 192 с.

37. Астапов В. Ф., Маркин А. А., Оленин С. И. Обоснование выбора мер конечных деформаций// Известия ТулГУ, сер. Математика, Механика, Информатика. -1995, т. 1, вып. 2, с. 7-11.

38. Морев П. Г. Учёт контакта инструмент-заготовка в конечноэлементных моделях// Известия ТулГУ, сер. Механика деформируемого твёрдого тела и обработка металлов давлением. 2003, вып. 1, с. 51-56.

39. Морев 77. Г. Конечноэлементный анализ упругопластических моделей в процессах ОМД// Сборник научных трудов ОрёлГТУ. 1996, т. 9, с. 42-47.

40. Морев 77. Г. Использование схемы Рунге-Кутта в конечноэлементном анализе процессов ОМД// Исследования в области теории, технологии и оборудования обработки металлов давлением. Межвузовский сб. научных трудов. Орёл: ОрёлГТУ-ТулГУ, 1998, с. 134-139

41. Морев 77. Г. Вариант метода конечных элементов для контактных задач с трением// Известия РАН, сер. Механика твёрдого тела. 2007, №4, с. 168-182

42. Ректорис К. Вариационные методы в математической физике и технике. М.: Мир, 1985,589 с.

43. Тимошенко С. 77., Гудъер Дж. Теория упругости. М.: Наука, 1979, 559 с.

44. Бураго Н. Г., Кукуджанов В. 77. Обзор контактных алгоритмов// Известия РАН, сер. Механика твёрдого тела. 2005, №1, с. 45-87

45. F. F. Mahmoud, S. S. Ali-Eldin, М. М. Hassan, S. А. Emam An incremental mathematical programming model for solving multi-phase contact problems// Computers and Structures, 1998, v. 68, p. 567-581

46. Wenhua Ling, H. K. Stolarski A contact algorithm for problems involving quadrilateral approximation of surfaces// Computers and Structures, 1997, v. 63,p. 963-975

47. G. Pietrzak, A. Curnier Large deformation frictional contact mechanics: continuum formulation and augmented Lagrangian treatment// Computer methods in applied mechanics and engineering, 1999, v. 177, p.351-381

48. Толоконников О. П., Маркин А. А., Астапов В. Ф. Исследование процесса формоизменения с учётом конечности деформаций// Прикладная механика, 1983, т.19, №10, с. 122-125

49. Бровко Г. Л. Материальные и пространственны^представления определяющих соотношений деформируемых сред// Прикладная математика и механика, 1990. т. 54, вып. 5. с. 814-824