3. Сафаров Д.Х. Эллиптическая регуляризация системы уравнений составного типа. - ДАН СССР, 1990, т. 311, №1, с. 36-39.
4. Сафаров Д.Х. Неклассические системы уравнений. - Душанбе: Дониш, 2008, 431 с.
5. Математическая энциклопедия. Полигармоническая функция. Т. 4. - М.: Советская энциклопедия, 1984, с. 403-404.
6. Векуа И.Н. Новые методы решения эллиптических уравнений. - М.-Л.: Гостехиздат, 1948, 296 с.
7. Миранда К. Уравнения с частными производными эллиптического типа. - М.: ИЛ, 1961, 256 с.
© Д.Х.Сафаров, С.С.Мирзоев, 2015
УДК: 539.37
Станкевич Игорь Васильевич
д.т.н., профессор МГТУ им. Н.Э. Баумана,
г. Москва, РФ E-mail: aplmex@yandex.ru
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ЗАДАЧ ТЕОРИИ УПРУГОСТИ С
ОДНОСТОРОННИМИ СВЯЗЯМИ
Аннотация
В работе рассматривается алгоритм построения численного решения контактной задачи теории упругости применительно к телу, которое имеет одностороннее контактное взаимодействие с абсолютно упругим полупространством.
Ключевые слова
Одностороннее контактное взаимодействие; контактная задача теории упругости; метод конечных
элементов.
Многие ответственные детали и узлы машиностроительных конструкций имеют контакт в пределах некоторой заданной поверхности [1]. В данной работе рассматривается частный случай контактного взаимодействия, когда нагруженное внешними силами упругое тело конечных размеров опирается на абсолютно жесткое полупространство (рис. 1, а). Контакт происходит по выделенной контактной поверхности, которая в общем случае может менять свои размеры. Процесс решения имеет итерационный характер и реализуются следующим образом. Исходная задач рассматривается как задача теории упругости,
при этом на контактной поверхности Sk задаются кинематические (по нормали П) и силовые (по
касательной Т) условия. Если силовые условия на Sk отсутствуют, то данная задача становится
стандартной задачей теории упругости. Силовые условия необходимо учитывать в том случае, когда задача решается с учетом трения [3-7]. Так как рассматривается односторонний контакт, то кинематические
условия в направлении нормали заранее известны, например, и\ = 0 (рис. 1), а силовые условия на
контактной поверхности Sk в касательном направлении после выполнения очередной итерации корректируются, для обеспечения выполнения принятого закона трения, и учитываются при проведении следующей итерации.
Для численного решения был использован метод конечных элементов (МКЭ) [2]. После минимизации
функционала полной потенциальной энергии линейно упругого тела, размещённого в пространстве П " и нагруженного массовыми и поверхностными, получается матричное уравнение (глобальная система линейных алгебраических уравнений), имеющее вид [2],
[К № } = {Щ. (1)
здесь [К] - глобальная матрица жесткости; {К} - глобальный вектор нагрузки, {и| - глобальный
вектор перемещений, состоящий из перемещений всех узлов конечно -элементной модели. При решении системы (1) не обходимо учесть заданные кинематические граничные условия, а также дополнительные
кинематические условия на поверхности , например, закрепление по нормали. В данной работе был использован метод сопряженных градиентов, что позволило учесть указанные условия непосредственно при реализации процедур метода сопряженных градиентов, не перестраивая для этого матрицу жесткости [ К ]
и вектор нагрузки {К} . Для учета касательных сил, возникающих на поверхности контакта из-за трения,
необходимо на каждой итерации вносить коррекцию в глобальный вектор нагрузки {К| [3, 4]. Алгоритм
коррекции касательных сил реализуется следующим образом. На первом шаге выполняется решение обычной упругой задачи без учета трения, то есть все узлы конечных элементов, расположенные на
контактной поверхности имели возможность при деформировании свободно перемещаться вдоль ограничивающей поверхности полупространства.
/ / / / / / / / / / / / / / / / / / / /
/ / / / / / / / У / / / / / / / / / / / / /
/ / / / / / / / / / / / / / / / / / / /
/ / / / / / / / / / / у / У / / / / / / / / /
/ / / / / / / / / / / / / / / / / /
/ / / / / / / /
/ / / / / / / / / / / / /
/ / / / / / / / / / / / / / / / / /
/ / / / / / / / / / / / / / / / /
/ \ / к X / к /
б
Рисунок 1 - Схема одностороннего контактного взаимодействия: а - заданное тело и кинематические связи; б - конечно-элементная модель
т т
Пусть и и и соответственно являются касательной и нормальной компонентами вектора
перемещения ит некоторого контактного узла Ат ^т = 1,М^ принадлежащего поверхности контакта »; здесь М -
конечно-элементной модели, число узлов конечно-элементной модели,
а
о т
принадлежащих поверхности контакта (рис.1). После вычисления касательного перемещения ЫТ ,
_т „ _т „ — т А
нормальной (Г и касательной (Г компонент вектора напряжении и в контактных узлах Ат
проверялось условие прилипания. Если устанавливалось, что
а
<М
ст„
то на следующей итерации в
рассматриваемом контактном узле Ат фиксировалось значение касательной компоненты Ы^ = 0 вектора перемещения ит, т.е. было запрещено горизонтальное перемещение рассматриваемого контактного узла
Ат . А в том случае, когда условие прилипания не выполнялось, проводилась коррекция касательной
рт
компоненты узловой силы г с помощью соотношения
pm (k+1) _ pm (k)
_mR(k) _mT (k)
т т
тах
а
nR(k)
а
где
pm
< k)
значение касательной компоненты узловой силы
г (к)
nT (k)
P П
-\Pm (k) ) '
(2)
в контактном узле
A
после
1 " Л Т>т (к) „ „ пп
выполнения к — ой итерации; Ар - заданное приращение касательной компоненты узловой силы Р
в контактном узле Am , например, APm (k)
а
PI
<k)
, где 0 < а < 1.
а
nR(k)
компоненты вектора напряжений, вычисленной по результатам к — ой итерации,
а
модуль касательной
чТ ( k )
теоретическое
а
nT (k)
м
а
(k)
а.
,(k)
значение модуля касательной компоненты вектора напряжений, здесь модуль нормальной компоненты вектора напряжений, вычисленной по результатам к — ой итерации.
т^т (к)
Вычисленное значение касательной компоненты узловой силы р вносится в соответствующую ячейку глобального вектора нагрузки {К} .
В выражении (2) при выборе знака необходимо руководствоваться следующим: если Г ( ) > 0 и
а
mR(k)
>
а
чТ (k)
, то выбирается знак «-»; если Г ( ) > 0 и
а
nT (k)
ö<
а
mR(k)
<
а
mT (k)
то
выбирается знак «+»; если P ( ) < 0
а
mR( k )
>
а
nT (k)
то выбирается знак «+»; если Г ( ) <0 и
а
T( k )
а
mR( k )
<
а
T( k )
то выбирается знак «-».
Так как численно добиться точного выполнения равенства, учитывающего скольжения при использовании закона Амонтона-Кулона, не удается, то вводится некоторый параметр 8 , который позволяет
организовать устойчивый численный итерационный процесс коррекции значений касательной компоненты пт „ т К
узловой силы р ив том случае, когда расчетные значения касательных напряжений о в
рассматриваемом контактном узле приближаются (по модулю) к теоретическим значениям 7^ снизу. В процессе выполнения итераций модуль разности значений касательных напряжений
Д т —тК(к) —тТ (к)
А О = 7Т — (7Т стремится к нулю.
Процесс вычислений заканчивается после оценки сходимости результатов итерации, например, для этого может быть использовано соотношение
и
s\l k+1 ^ = max
m=1,M
m (k+1) m (k)
u u
<e„> (3)
где и " (к+1 - вектор перемещения ит контактного узла Ат , после выполнения (к + 1) — ой
итерации; £ - требуемая точность по перемещениям контактных узлов.
В качестве примера рассмотрено плоское напряженное состояние пластинки, занимающей в
двухмерном пространстве 0 с декартовой системой координат О X | X 0 область О с размерами в плане 0,11м X 0,05 м. Пластинка выполнена из материала с модулем упругости Е = 200 ГПа и коэффициентом Пуассона 0,3. На нижней поверхности » имеются выступы треугольной формы (рис. 1). Внешними
вершинами выступов Ат , т = 1, М, где М - число вершин выступов (М = 5), пластинка опирается на абсолютно жесткое полупространство. Роль контактной поверхности » выполняет поверхность полупространства. Поверхность » , на которой заданы кинематические граничные условия (2), состоит из трех участков , »12 и : на участке заданы отрицательные перемещения и2 = - 0,00004м в направлении оси Ох2 , а на участке »12 положительные перемещения и, = + 0,0001м в направлении оси Ох,. Поверхность »13 зафиксирована от горизонтальных перемещений (и, — 0), но вертикальные перемещения и2 точек поверхности допустимы. На поверхности заданы условия не проникновения в
полупространство, то есть принято, что для вершин выступов Ат и^ = 0 , т = 1, М, а ограничений на
„. т
горизонтальные перемещения нет.
На рис. 2 представлены графики изменения компоненты С22 вектора напряжений Гт и касательной
компоненты и1 вектора перемещения ит контактных узлов Ат (т = 1,5) в зависимости от числа
итераций N. Из графиков видно, что использование соотношения (2) для коррекции касательных сил, учитывающих эффект трения в контактных узлах конечно-элементной модели, позволяет получить устойчивое сходящееся решение задачи теории упругости с односторонним контактом.
б
Рисунок 3 - Изменения значений компоненты С22 вектора напряжений Гт (а) и касательной
^^ т а
компоненты и, = ЫТ вектора перемещения и (б) контактных узлов (см. рис. 1, б) в зависимости от
числа итераций N
а