A nonlinear boundary integral equations method for the solving of quasistatic elastic contact problem with Coulomb friction

Full Text

Введение. Для решения контактных задач теории упругости с учетом трения часто используются численные методы, основанные на вариационной постановке задачи [1, 2] либо на её сведении к различным нелинейным операторным уравнениям [3-7]. Основная трудность реализации вариационных методов заключается в необходимости рассматривать сложные задачи нелинейного программирования. Иногда эту трудность удается устранить путем использования нелинейных операторных уравнений для моделирования контактного взаимодействия тел. Однако в большинстве случаев при использовании таких уравнений авторы либо не учитывают трение [3-5], либо его учет осуществляется при упрощенных граничных условиях, соответствующих полному проскальзыванию тел [6, 7]. Хотя метод нелинейных граничных интегральных уравнений, предложенный в работе [8], и позволяет рассматривать закон трения Кулона в сложном классическом виде, но при использовании этого метода трудно выразить условия единственности решения контактной задачи и доказать сходимость итерационных процессов, применяемых для получения этого решения [9]. В настоящей статье предложена модификация нелинейных граничных интегральных уравнений [8], позволяющая избавиться от указанных трудностей. Данная модификация основана на упрощении граничных условий квазистатической контактной задачи на каждом шаге процесса нагружения тел. Это упрощение позволяет свести контактную задачу к последовательному решению серии однотипных нелинейных операторных уравнений второго рода с квазисжимающим оператором в правой части. Такая особенность этих уравнений позволяет разработать эффективный метод их приближённого решения. Целью данной статьи является разработка такого метода и получение строгого математического обоснования основных его этапов. 1. Постановка контактной задачи. Предположим, что каждое из двух тел на некоторой части своей поверхности сцеплено с абсолютно твёрдым телом, которое будем называть жесткой опорой. Первоначально касаясь в одной точке, тела входят в контакт за счет того, что опора верхнего тела совершает заданное поступательное перемещение γ (t), зависящее от времени t, в то время как опора нижнего тела остается неподвижной. Требуется найти поверхность соприкасания тел после установившего равновесия и распределение контактной нагрузки на этой поверхности. Относительно взаимодействующих тел принимаются следующие допущения: 1) оба тела являются линейно-упругими и изотропными; 2) упругие перемещения точек малы по сравнению с размерами поверхности контакта; 3) в каждой точке поверхности каждого из тел в зоне возможного контакта существует единственная касательная плоскость к этой поверхности; 4) операторы влияния поверхностных нагрузок на поверхностные упругие перемещения для каждого из тел являются известными и допускают интегральное представление. Будем считать, что вектор-функция γ (t) медленно изменяется с возрастанием t от 0 до фиксированного значения T > 0, и пренебрегать инерционными и волновыми эффектами. Предположим, что процесс взаимодействия тел сопровождается поверхностным трением, подчиняющемся закону Кулона [10]. Будем ещё полагать, что в любой момент времени t > 0 поверхность соприкасания тел является плоской и содержится в известной ограниченной области Ω, принадлежащей общей для тел касательной плоскости Π, проходящей через точку их начального касания. Граница областей сцепления и проскальзывания заранее не известна и подлежит определению в ходе решения. Введем прямоугольную систему координат Ozxy c началом в точке первоначального касания тел и осью z, направленной внутрь верхнего тела перпендикулярно плоскости Π. Будем рассматривать процесс нагружения тел лишь в фиксированные моменты времени t0 , t1 , t2 , ..., tl такие, что 0 = t0 < t1 < t2 < · · · < tl = T . Будем также предполагать, что при каждом значении t = ti для системы контактирующих тел устанавливается равновесие. Таким образом, процесс нагружения тел происходит пошагово и состоит из конечного числа l последовательных состояний равновесия. Пусть индекс i определяет номер шага нагружения и вектор-функция γ (t) принимает на каждом шаге нагружения значение γi = γ (ti ) = (∆1i , ∆2i , ∆3i ) (рис. 1). Граничные условия контактного взаимодействия тел на i-том шаге нагружения будем выражать с помощью двух вектор-функций pi (s) = p1i (s), p2i (s), p3i (s) , vi (s) = v1i (s), v2i (s), v3i (s) , где s = (x, y) - произвольная точка области Ω, определяющая пару противолежащих точек (рис. 1) на поверхностях взаимодействующих тел (предполагается, что произвольная точка поверхности верхнего тела может вступить в контакт только с противолежащей точкой поверхности нижнего тела). Вектор-функция pi (s) определяет распределение удельной контактной нагрузки, передаваемой от нижнего тела к верхнему телу через точки области Ω на i-том шаге нагружения. Вектор-функция vi (s) задает перемещения точек поверхности верхнего тела относительно противолежащих точек поверхности нижнего тела на i-том шаге нагружения. Зависимость компонент этих вектор-функций от времени выражается лишь в дискретном изменении их значений при переходе к следующему шагу нагружения. Отметим, что значения 1, 2, 3 индексов соответствуют направлениям координатных осей Oz, Ox, Oy. С учетом принятых предположений граничные условия контактного взаимодействия тел можно выразить следующей системой неравенств и уравнений: p1i (s) 0; p22i (s) + v1i (s) 2 p3i (s) 0; p1i (s)v1i (s) = 0; µp1i (s); [v˙ 2i (s)]2 + [v˙ 3i (s)]2 p2i (s) + µp1i (s)v˙ 2i (s) = 0; (1) [v˙ 2i (s)]2 + [v˙ 3i (s)]2 p3i (s) + µp1i (s)v˙ 3i (s) = 0, где s ∈ Ω, i = 1, l; положительная константа µ - коэффициент трения, а v˙ 2i (s), v˙ 3i (s) - компоненты вектора скорости относительного проскальзывания тел в точке s ∈ Ω на i-том шаге нагружения. Первое из соотношений (1) выражает знакопостоянство контактного давления, второе означает отсутствие взаимного проникания тел, третье означает отсутствие контактного давления за пределами площадки контакта. Последние три соотношения системы (1) выражают закон трения Кулона. Будем считать, что компоненты pki , vki вектор-функций pi (s), vi (s) принадлежат гильбертову пространству L2 (Ω), а сами вектор-функции pi (s), vi (s) являются элементами гильбертова пространства L32 (Ω) [9]. Связь между вектор-функциями pi (s) и vi (s) выражается соотношениями 3 Akj (pji )s - fki (s) ≡ Fk (pi , ∆ki )s ∀k = 1, 3, vki (s) = j=1 (2) Рис. 1. Схема контактного взаимодействия упругих тел [Figure 1. The scheme of contact interaction of elastic bodies] где Akj : L2 (Ω) → L2 (Ω) - операторы влияния поверхностных нагрузок на поверхностные упругие перемещения; fi (s) = (f1i (s), f2i (s), f3i (s)) ∈ L32 (Ω) - известная вектор-функция, для которой f1i (s) = -δ0 (s) - ∆1i , f2i (s) = -∆2i , f3i (s) = -∆3i , а функция δ0 (s) ∈ L2 (Ω) задает первоначальный зазор между телами, определяемый уравнениями их поверхностей. Выразим компоненты v˙ 2i (s), v˙ 3i (s) вектора скорости относительного проскальзывания тел в точке s ∈ Ω на i-том шаге нагружения с помощью конечных разностей: v˙ ki (s) = vk i (s) - vk i-1 (s) 1 ˜ ki )s , = Fk (pi , ∆ hi hi i = 1, l, k = 2, 3, (3) ˜ ki (s) = ∆ki - Fk (pi-1 , ∆ki-1 ) , а Fk (pi , ∆ki ) определяется где hi = ti - ti-1 , ∆ s s соотношениями (2). Тогда систему (1) можно записать в виде p1i (s) 0; F1 (pi , ∆1i )s p22i (s) + p23i (s) 0; p1i (s)F1 (pi , ∆1i )s = 0; µp1i (s); ˜ 2i )s + F 2 (pi , ∆ ˜ 3i )s p2i (s) + µp1i (s)F2 (pi , ∆ ˜ 2i )s = 0; F22 (pi , ∆ 3 (4) ˜ 2i )s + F 2 (pi , ∆ ˜ 3i )s p3i (s) + µp1i (s)F3 (pi , ∆ ˜ 3i )s = 0; F22 (pi , ∆ 3 где s ∈ Ω, i = 1, l. Таким образом, рассматриваемая контактная задача на i-том шаге нагружения состоит в определении неизвестных функций pki (s) ∈ L2 (Ω), k = 1, 3, которые почти всюду на Ω удовлетворяют соотношениям (4). Отметим, что из соотношений (3) следует, что для нахождения решения системы (4) на i-том шаге нагружения должно быть известно решение этой системы на (i - 1)-м шаге. В качестве начальных условий на нулевом шаге (в момент начального касания тел) примем тождественное равенство нулю всех компонент векторфункций γ0 и p0 (s). Для упрощения постановки контактной задачи рассмотрим следующую модификацию соотношений, выражающих закон трения Кулона в системе (4). Вначале запишем эту систему в эквивалентном виде: p1i (s) 0; p22i (s) + F1 (pi , ∆1i )s 0; 2 p3i (s) µ (p1i-1 (s) p1i (s)F1 (pi , ∆1i )s = 0; + εi (s)) ; ˜ 2i )s + F 2 (pi , ∆ ˜ 3i )s p2i (s) + µ (p1i-1 (s) + εi (s)) F2 (pi , ∆ ˜ 2i )s = 0; F22 (pi , ∆ 3 (5) ˜ 2i )s + F 2 (pi , ∆ ˜ 3i )s p3i (s) + µ (p1i-1 (s) + εi (s)) F3 (pi , ∆ ˜ 3i )s = 0, F22 (pi , ∆ 3 где s ∈ Ω, i = 1, l; εi (s) = p1i (s) - p1i-1 (s). Затем будем предполагать, что при достаточно большом числе шагов нагружения l разности p1i (s)-p1i-1 (s) малы по сравнению с величинами контактных давлений p1i (s) на каждом i-том шаге нагружения. Наконец, учитывая последнее предположение, отбросим в системе (5) все приращения εi (s). Тогда система (5) примет следующий модифицированный вид: p1i (s) 0; F1 (pi , ∆1i )s p22i (s) + p23i (s) 0; p1i (s)F1 (pi , ∆1i )s = 0; µp1i-1 (s); ˜ 2i )s + F 2 (pi , ∆ ˜ 3i )s p2i (s) + µp1i-1 (s)F2 (pi , ∆ ˜ 2i )s = 0; F22 (pi , ∆ 3 (6) ˜ 2i )s + F 2 (pi , ∆ ˜ 3i )s p3i (s) + µp1i-1 (s)F3 (pi , ∆ ˜ 3i )s = 0, F22 (pi , ∆ 3 где s ∈ Ω, i = 1, l. При использовании указанной модификации закона трения Кулона рассматриваемая контактная задача на i-том шаге нагружения сводится к определению неизвестных функций p1i (s), p2i (s), p3i (s) ∈ L2 (Ω), которые почти всюду на Ω удовлетворяют соотношениям (6). Принятую модификацию закона трения Кулона, выраженную последними тремя соотношениями системы (6), можно истолковать как наличие запаздывания в выражениях µp1i (s), входящих в последние три соотношения системы (6) и ограничивающих модуль удельной касательной контактной нагрузки. При этом очевидно, что на i-том шаге нагружения ограничивающее контактное давление p1i-1 (s) является известным, что существенно упрощает граничные условия контактной задачи. Если система (6) решена для i-того шага нагружения, то конфигурацию области контакта на этом шаге можно определить следующим образом: (i) Ω0 = s ∈ Ω | F1 (pi , ∆1i )s = 0 , (i) (i) ˜ 2i )s + F 2 (pi , ∆ ˜ 3i )s = 0 , ΩC = s ∈ Ω0 | F22 (pi , ∆ 3 (i) (i) ˜ 2i )s + F 2 (pi , ∆ ˜ 3i )s > 0 , ΩΠ = s ∈ Ω0 | F22 (pi , ∆ 3 (i) (i) (i) где Ω0 - область контакта, ΩC - зона сцепления, ΩΠ - зона проскальзывания. Отметим, что при изменении индекса i от 1 до l соотношения (6) представляют собой серию систем, подлежащих последовательному решению. Окончательным решением рассматриваемой квазистатической задачи усло(l) (l) (l) вимся считать функции pkl (s), k = 1, 3 и области Ω0 , ΩC , ΩΠ , найденные на последнем шаге нагружения l. 2. Операторное уравнение контактной задачи. Для неизвестных векторфункций p1 (s), p2 (s), . . . , pl (s) ∈ L32 (Ω) серия систем соотношений (6) эквивалентна серии нелинейных операторных уравнений [9] pi = Gµp1 i-1 pi - E A (pi ) - f˜i , i = 1, l, (7) в которых E - произвольная положительная константа, f˜i ∈ L32 (Ω) и определяются равенствами ˜ 2i (s), -∆ ˜ 3i (s) . f˜i = -δ0 (s) - ∆1i , -∆ Линейный ограниченный оператор влияния A : L32 (Ω) → L32 (Ω), входящий в уравнения (7), задается соотношениями x = (x1 , x2 , x3 ) , y = (y1 , y2 , y3 ) ∈ L32 (Ω); 3 (8) Akj (xj ) ∀k = 1, 3. y = A(x), yk = j=1 Входящий в правую часть соотношений (7) нелинейный непрерывный оператор Gµg : L32 (Ω) → L32 (Ω) определяется равенствами x = (x1 , x2 , x3 ), y = (y1 , y2 , y3 ) ∈ L32 (Ω); y = Gµg (x), y1 (s) = h x1 (s) , y2 (s) = q x2 (s), x3 (s), µg(s) , y3 (s) = q x3 (s), x2 (s), µg(s) , s ∈ Ω, (9) в которых g(s) ∈ L2 (Ω) - заданная, почти всюду на Ω неотрицательная функция, а функции h(x) и q(x, y, z) определяются соотношениями x, если x 0; 0, если x < 0;   x, если  q(x, y, z) = xz  , если  x2 + y 2 h(x) = x2 + y 2 z; (10) x2 + y 2 > z. Для получения оператора Gµp1 i-1 : L32 (Ω) → L32 (Ω), который фактически используется в выражении правой части уравнения (7), нужно в соотношениях (9), задающих оператор Gµg : L32 (Ω) → L32 (Ω), положить g(s) = p1 i-1 (s) почти всюду на Ω. Важной особенностью [9] каждого из уравнений (7) является то, что множество его решений в L32 (Ω) не изменяется при изменении входящего в это уравнение положительного параметра E. Оператор Gµg : L32 (Ω) → L32 (Ω), определяемый соотношениями (9), является метрическим проектором [11] гильбертова пространства L32 (Ω) на выпуклое замкнутое множество U = (x1 , x2 , x3 ) ∈ L32 (Ω) | x1 (s) 0, x22 (s) + x23 (s) µg(s), s ∈ Ω этого пространства. Этот факт следует из очевидных свойств функций h(x) и q (x, y, z), определяемых соотношениями (10). Действительно, функция h(x) является метрическим проектором одномерного евклидова пространства R1 на промежуток [0; +∞) ⊂ R1 , а отображение ϕr : R2 → R2 вида x = (x1 , x2 ), y = (y1 , y2 ) ∈ R2 ; y = ϕr (x); y1 = q(x1 , x2 , r); y2 = q(x2 , x1 , r); r 0, осуществляет метрическую проекцию двумерного евклидова пространства R2 на замкнутый круг (x1 , x2 ) ∈ R2 | x21 + x22 r этого пространства. Отметим, что входящая в (7) известная неотрицательная функция p1 i-1 ∈ L2 (Ω) определяется из решения pi-1 ∈ L32 (Ω) этого уравнения на (i - 1)-м шаге нагружения, а известный элемент f˜i ∈ L32 (Ω) также зависит от решения pi-1 . Таким образом, решение pi ∈ L32 (Ω) уравнения (7) при фиксированном значении индекса i можно получить после нахождения решений p1 , p2 , . . . , pi-1 ∈ L32 (Ω) этого уравнения на предшествующих шагах нагружения (при получении p1 предполагается, что в момент начального касания тел вектор-функция контактной нагрузки и жесткие сближения тел тождественно равны нулю на Ω). Окончательным решением квазистатической контактной задачи, описываемой серией уравнений (7), будем считать решение pl ∈ L32 (Ω), полученное на последнем l-том шаге нагружения. Для многих важных для практики случаев входящие в соотношения (8) линейные ограниченные операторы Akj : L2 (Ω) → L2 (Ω) допускают интегральное представление. В этих случаях каждое из операторных уравнений (7) при фиксированном значении индекса i представляет собой систему трех нелинейных интегральных уравнений относительно трех неизвестных функций p1 i (s), p2 i (s), p3 i (s), задающих распределение удельной контактной нагрузки, передаваемой от одного тела к другому. 3. Единственность решения контактной задачи. Для доказательства единственности решения контактной задачи, описываемой серией операторных уравнений (7), достаточно доказать единственность решения каждого из уравнений этой серии. Фиксируя индекс i в (7) и обозначая p = pi ∈ L32 (Ω), g = p1 i-1 ∈ L2 (Ω), b = f˜i ∈ L32 (Ω), получим более простой вид каждого операторного уравнения из совокупности (7): p = Gµg p - E A(p) - b . (11) В дальнейшем будем полагать, что линейный ограниченный оператор влияния A : L32 (Ω) → L32 (Ω), входящий в уравнение (11), является самосопряженным и положительным: A(x), y = x, A(y) ∀x, y ∈ L32 (Ω), A(x), x > 0 ∀x ∈ L32 (Ω)\{θL32 (Ω) }, 312 (12) где (x, y) - скалярное произведение элементов x, y пространства L32 (Ω), θL32 (Ω) - нулевой элемент пространства L32 (Ω). Самосопряженность оператора A(p) обусловлена свойством взаимности упругой среды, а его положительность гарантирована положительностью энергии упругих деформаций системы взаимодействующих тел, подверженных ненулевому нагружению. Сформулируем и докажем теорему о единственности решения уравнения (11) при условии, что оператор влияния A удовлетворяет соотношениям (12). Теорема 1. Пусть линейный ограниченный оператор A : L32 (Ω) → L32 (Ω) удовлетворяет условиям (12). Тогда для любого элемента b ∈ L32 (Ω) любой заданной и почти всюду неотрицательной на Ω функции g ∈ L2 (Ω) и любого действительного числа µ > 0 уравнение (11) не может иметь более одного решения в пространстве L32 (Ω). Д о к а з а т е л ь с т в о. Предположим противное. Пусть p и p˜ - два различных решения уравнения (11) в пространстве L32 (Ω). Для оценивания нормы разности элементов p и p˜ воспользуемся нерастягиваемостью на этом пространстве метрического проектора Gµg : L32 (Ω) → L32 (Ω) [11]: p - p˜ = Gµg p - E A(p) - b - Gµg p˜ - E A (˜ p) - b p - p˜ - EA (p - p˜) . Возведя обе части полученного неравенства в квадрат и перейдя в нём от норм элементов к скалярным произведениям, получим оценку p - p˜ 2 p - p˜ 2 - 2E p - p˜, A (p - p˜) + E 2 A (p - p˜) 2 , означающую справедливость неравенства p - p˜, A (p - p˜) E A (p - p˜) 2 2 . Очевидно, что полученное неравенство будет выполняться для любого положительного значения E. Переходя в этом неравенстве к пределу при E → +0, получим соотношение p - p˜, A p - p˜ 0, которое противоречит второму из условий в (12). Полученное противоречие означает, что доказываемое утверждение верно. Доказанная теорема означает, что при выполнении её условий рассматриваемая квазистатическая контактная задача, описываемая серией уравнений (7), не может иметь более одного решения в пространстве L32 (Ω). 4. Регуляризация операторного уравнения контактной задачи. Разработка численного алгоритма решения уравнения (11) осложняется тем, что во многих случаях входящий в него линейный оператор влияния A : L32 (Ω) → L32 (Ω) является вполне непрерывным [9, 11]. Это обстоятельство приводит к тому, что задача нахождения неизвестного элемента p ∈ L32 (Ω) по известным элементам b ∈ L32 (Ω), g ∈ U0 из решения уравнения (11) является некорректно поставленной в смысле А. Н. Тихонова [12] (здесь и ниже символом U0 обозначен конус неотрицательных функций пространства L2 (Ω) [13]). Поэтому 313 С т р е л я е в Ю. М. вместо уравнения (11) целесообразно рассматривать его регуляризированный аналог p = Gµg p - E εp + A(p) - b , (13) в котором ε > 0 есть некоторое фиксированное число - параметр регуляризации. Сформулируем и докажем две теоремы, из которых следует, что уравнение (13) является регуляризированным аналогом уравнения (11). Теорема 2. Пусть линейный ограниченный оператор A : L32 (Ω) → L32 (Ω) удовлетворяет условиям (12). Тогда для любого фиксированного положительного числа ε задача нахождения неизвестного элемента p ∈ L32 (Ω) по известным элементам b ∈ L32 (Ω), g ∈ U0 из решения уравнения (13) является корректно поставленной в смысле А. Н. Тихонова. Д о к а з а т е л ь с т в о. Для доказательства теоремы достаточно показать, что для любых элементов b ∈ L32 (Ω), g ∈ U0 уравнение (13) однозначно разрешимо в пространстве L32 (Ω) и его решение p ∈ L32 (Ω) непрерывно зависит от элементов b ∈ L32 (Ω), g ∈ U0 [12]. Пусть ε - фиксированное положительное число и b ∈ L32 (Ω), g ∈ U0 . Выберем положительное число E в уравнении (13) так, чтобы выполнялось неравенство E< 1 ε+ A (14) , ∗ где A ∗ - норма оператора A : L32 (Ω) → L32 (Ω) [14]. Докажем, что из неравенства (14) следует оценка I - E (εI + A) ∗ (15) 1 - Eε < 1, где I - тожественный оператор, действующий в пространстве L32 (Ω). Из неравенства (14) и условий (12) следует, что линейный ограниченный оператор (I - E (εI + A)) : L32 (Ω) → L32 (Ω) является самосопряженным и положительным. Найдя его норму с помощью известного соотношения [14], получим для этой нормы очевидную оценку: I - E (εI + A) ∗ = sup x - E εx + A(x) , x = x =1 = 1 - Eε - E inf A(x), x 1 - Eε = q < 1, x =1 из которой следует неравенство (15). Из этого неравенства и нерастягиваемости метрического проектора Gµg : L32 (Ω) → L32 (Ω) следует сжимаемость оператора F : L32 (Ω) → L32 (Ω), стоящего в правой части уравнения (13). Действительно, для любых элементов x, y ∈ L32 (Ω) можно получить оценки F (x) - F (y) = = Gµg x - E εx + A(x) - b - Gµg y - E εy + A(y) - b x - y - E ε (x - y) + A (x - y) I - E (εI + A) ∗ · x - y 314 q x-y , Метод нелинейных граничных интегральных уравнений. . . где 0 < q < 1. Таким образом, однозначная разрешимость уравнения (13) в пространстве L32 (Ω) для любых элементов b ∈ L32 (Ω), g ∈ U0 следует из принципа сжимающих отображений [15]. Пусть элементы p , p ∈ L32 (Ω) являются решениями уравнения (13), полученными при b = b ∈ L32 (Ω), g = g ∈ U0 и b = b ∈ L32 (Ω), g = g ∈ U0 соответственно. Тогда, учитывая нерастягиваемость оператора Gµg : L32 (Ω) → L32 (Ω) и соотношения (10), (15), можно получить неравенства p -p = = Gµg p - E εp + A(p ) - b - Gµg = Gµg p - E εp + A(p ) - b +Gµg p - E εp + A(p ) - b Gµg p - E εp + A(p ) - b p - E εp + A(p ) - b - Gµg p - E εp + A(p ) - b - Gµg p - E εp + A(p ) - b - Gµg · p -p +E b -b (1 - Eε) p - p из которых для нормы p - p p -p + p - E εp + A(p ) - b (p - p ) - E ε(p - p ) + A(p - p ) + E(b - b ) + µ g - g ∗ + - Gµg p - E εp + A(p ) - b + Gµg p - E εp + A(p ) - b I - E (εI + A) = +µ g -g +E b -b 1 1 +µ g -g 1, следует оценка 1 b -b ε + µ g -g Eε 1. Здесь и ниже символами · и · 1 обозначены нормы элементов пространств L32 (Ω) и L2 (Ω) соответственно. Эта оценка, полученная при фиксированном ε > 0 для любых элементов b , b ∈ L32 (Ω) и g , g ∈ U0 , означает, что решение уравнения (13) в пространстве L32 (Ω) непрерывно зависит от элементов b ∈ L32 (Ω), g ∈ U0 . Теорема 3. Пусть линейный вполне непрерывный оператор A : L32 (Ω) → 3 L2 (Ω) удовлетворяет условиям (12) и для некоторых фиксированных элементов b ∈ L32 (Ω), g ∈ U0 уравнение (11) имеет решение p∗ ∈ L32 (Ω). Тогда lim ε→+0 p∗ - pε = 0, где pε ∈ L32 (Ω) является решением уравнения (13), полученным при заданном положительном значении параметра регуляризации ε. Д о к а з а т е л ь с т в о. Из очевидных равенств p∗ = Gµg (p∗ - E (A(p∗ ) - b)) , pε = Gµg pε - E εpε + A pε - b вытекает оценка p∗ - pε 2 p∗ - pε - E (A(p∗ - pε ) - εpε ) 2 , справедливая для любого числа E > 0. Переходя в этой оценке от норм элементов к скалярным произведениям, получим неравенство p∗ - pε , A (p∗ - pε ) - εpε E A (p∗ - pε ) - εpε 2 2 . 315 С т р е л я е в Ю. М. Переходя в этом неравенстве к пределу при E → +0, получим условие p∗ - pε , A (p∗ - pε ) - εpε (16) 0. Используя положительность оператора A и неравенство Коши-Буняковского [15], получим из условия (16) оценки 0 p∗ - pε , A (p∗ - pε ) (p∗ - pε , εpε ) ε pε p∗ - pε , влекущие за собой неравенство pε p∗ . (17) Очевидно, что неравенства (16) и (17) выполняются для любого числа ε > 0. Предположим теперь, что {εn } - произвольная последовательность положительных чисел, для которой lim εn = 0. n→∞ Введя обозначение pεn = pn , покажем, что lim pn - p∗ = 0. n→∞ Это и будет означать справедливость доказываемой теоремы. Для доказательства сильной сходимости последовательности {pn } элементов L32 (Ω) к элементу p∗ ∈ L32 (Ω) достаточно показать, что из любой подпоследовательности {pnk } последовательности {pn } можно извлечь другую подпоследовательность {pnkj }, которая сильно сходится к элементу p∗ . Покажем это. Пусть {pnk } - произвольная подпоследовательность последовательности {pn }. Поскольку подпоследовательность {pnk } является ограниченной в силу неравенства (17), существует такая её подпоследовательность {pnkj } и такой элемент p˜ ∈ L32 (Ω), что {pnkj } слабо сходится в L32 (Ω) к p˜ (см. [11, теорема 1.7]). Из оценки (16) при ε = εnkj вытекает неравенство p∗ - pnkj , A(p∗ - pnkj ) εnkj (p∗ , pnkj ) - εnkj pnkj 2 , справедливое для всех значений индекса nkj . Переходя в этом неравенстве к пределу при nkj → ∞, получим неравенство p∗ - p˜, A (p∗ - p˜) 0, которое в силу второго из условий (12) означает, что p∗ = p˜. Следовательно, {pnkj } слабо сходится к элементу p∗ . Но поскольку ещё из неравенства (17) следуют оценки pnkj p∗ ∀nkj , подпоследовательность {pnkj } сильно сходится к p∗ (см. [11, теорема 1.4]). 316 Метод нелинейных граничных интегральных уравнений. . . Из доказанных теорем 2 и 3 следует, что уравнение (13) является регуляризирующим аналогом уравнения (11). 5. Дискретизация регуляризированного уравнения. Наличие слагаемого εp в уравнении (13) позволяет построить эффективный устойчивый численный алгоритм приближенного решения этого уравнения. Такой алгоритм основывается на следующей аппроксимационной теореме. Пусть b ∈ L32 (Ω), g ∈ U0 , а линейный ограниченный опе3 ратор A : L2 (Ω) → L32 (Ω) удовлетворяет условиям (12). Пусть последовательность {An } линейных ограниченных операторов, действующих в L32 (Ω), а также последовательность {bn } элементов L32 (Ω) и последовательность {gn } элементов U0 удовлетворяют следующим условиям: Теорема 4. lim A - An ∗ = 0, lim b - bn = 0, n→∞ lim g - gn 1 = 0. n→∞ (18) n→∞ Тогда для любого числа ε > 0 существует натуральное число n0 , зависящее только от ε, такое, что для всех натуральных n n0 уравнение (19) p = Gµgn p - E εp + An (p) - bn является однозначно разрешимым в L32 (Ω) и справедливо равенство lim pn - pε = 0, n→∞ где pn , pε ∈ L32 (Ω) являются решениями уравнений (19) и (13) соответственно. Д о к а з а т е л ь с т в о. Выберем для заданного числа ε > 0 номер n0 и число E > 0 в уравнении (19) так, чтобы выполнялись следующие неравенства: A - An ∗ < ε ∀n 2 n0 , E< 1 ε+ A . ∗ Используя эти неравенства, можно получить оценку I - E (εI + An ) 1-E ∗ ε < 1 ∀n 2 (20) n0 . Действительно, I - E (εI + An ) ∗ I - E (εI + A) = I - E (εI + An - A + A) ∗ + E A - An ∗ ∗ sup x - E εx + A(x) , x + E x =1 1 - Eε + E ε 2 ε ε = 1 - E < 1. 2 2 Из неравенства (20) следует сжимаемость оператора Fn : L32 (Ω) → L32 (Ω), стоящего в правой части уравнения (19) (см. доказательство теоремы 2), которая 317 С т р е л я е в Ю. М. влечет за собой однозначную разрешимость этого уравнения в пространстве L32 (Ω) для всех натуральных n n0 . Докажем теперь сильную сходимость в пространстве L32 (Ω) последовательности {pn } решений уравнений (19) к решению pε ∈ L32 (Ω) уравнения (13). Для этого рассмотрим последовательность {αn } элементов L32 (Ω), определяемую равенством αn = An (pn ) - A (pn ) + b - bn , n n0 . Учитывая условия (18) и ограниченность в L32 (Ω) последовательности {pn }, имеющую место в силу ограниченности в L32 (Ω) сходящейся последовательности {bn }, а также неравенства (20), получим для αn оценки αn An - A ∗ pn + b - bn An - A ∗ · β˜ + b - bn , где β˜ - некоторое положительное число. Из последних оценок следует, что (21) lim αn = 0. n→∞ Поскольку из теоремы 2 следует непрерывная зависимость решения уравнения (13) от элементов b ∈ L32 (Ω), g ∈ U0 , из очевидных соотношений pε = Gµg pε - E εpε + A (pε ) - b , pn = Gµgn pn - E εpn + A (pn ) - b + αn , а также из последнего равенства (18) и доказанного соотношения (21) следует, что lim pn - pε = 0. n→∞ Из теорем 1-4 следует, что приближенное решение рассматриваемой контактной задачи на каждом шаге нагружения можно получить, решив уравнение (19) при достаточно малом значении ε > 0 и достаточной близости оператора An и элементов bn , gn к оператору A и элементам b, g соответственно. Построим такие удовлетворяющие условиям теоремы 4, аппроксимирующие последовательности {An }, {bn } и {gn }, для которых уравнение (19) можно решить численно с любой степенью точности. Для этого зададим область Ω в виде открытого квадрата, ограниченного прямыми, параллельными осям некоторой декартовой системы координат, введенной на общей для взаимодействующих тел касательной плоскости. Далее для каждого натурального n разобьем область Ω на n2 непересекающихся квадратных областей ω1 , ω2 , . . . , ωn2 равной площади, ориентированных подобно квадрату Ω. Будем полагать, что в выражениях (8) линейные ограниченные операторы Aij : L2 (Ω) → L2 (Ω) допускают интегральное представление Kij (t, s)x(s)ds, t ∈ Ω, i, j = 1, 3 Aij (x)t = Ω и являются вполне непрерывными. Зададим линейный вполне непрерывный оператор An : L32 (Ω) → L32 (Ω) следующим образом: 318 Метод нелинейных граничных интегральных уравнений. . . x = (x1 , x2 , x3 ) , y = (y1 , y2 , y3 ) ∈ L32 (Ω); 3 (n) y = An (x), yi = (n) Aij (xj )t = Ω Aij (xj ) , i = 1, 3; j=1 (n) Kij (t, s)xj (s)ds, (22) t ∈ Ω, i, j = 1, 3. (n) Входящие в соотношения (22) функции Kij (t, s) определяются равенствами (n) Kij (t, s) = 1 mes(ωk ) mes(ωm ) Kij (t, s)ds dt, ωk ωm если t ∈ ωk , s ∈ ωm ; k, m = 1, n2 ; i, j = 1, 3, (23) в которых mes(ωk ) есть площадь квадрата ωk . Элементы bn ∈ L32 (Ω) и gn ∈ U0 зададим в виде bn = Pn (˜b1 ), Pn (˜b2 ), Pn (˜b3 ) , gn = Pn (g), где ˜b1 , ˜b2 , ˜b3 = b и оператор Pn : L2 (Ω) → L2 (Ω) определяется соотношениями Pn (x)t = 1 mes(ωk ) t ∈ ωk , k = 1, n2 . x(s)ds, ωk Очевидно, что в случае полной непрерывности линейных интегральных операторов Aij : L2 (Ω) → L2 (Ω) построенные аппроксимирующие последовательности {An }, {bn } и {gn } удовлетворяют условиям (18). Кроме этого, эти последовательности заданы так, что оператор Fn : L32 (Ω) → L32 (Ω), стоящий в правой части уравнения (19), отображает пространство L32 (Ω) в конечномерное подпространство Hn = Pn (L2 (Ω)) × Pn (L2 (Ω)) × Pn (L2 (Ω)) пространства L32 (Ω). Следовательно, каждое решение pn ∈ L32 (Ω) уравнения (19) содержится в Hn и может быть найдено из решения следующей системы 3n2 скалярных уравнений с 3n2 неизвестными: 3n2 x3k-2 = h x3k-2 - E a3k-2j xj - b3k-2 ; j=1 3n2 x3k-1 = q x3k-1 - E a3k-1j xj - b3k-1 , j=1 3n2 x3k - E a3kj xj - b3k , µgk ; j=1 (24) 3n2 a3kj xj - b3k , x3k = q x3k - E j=1 3n2 x3k-1 - E a3k-1j xj - b3k-1 , µgk ; j=1 k = 1, n2 . 319 С т р е л я е в Ю. М. Решение pn = (p1n , p2n , p3n ) ∈ Hn уравнения (19) выражается через решение x1 , x2 , . . . , x3n2 системы (24) соотношениями x3k-2 = p1n (s), x3k-1 = p2n (s), x3k = p3n (s) при s ∈ ωk ∀k = 1, n2 . Неотрицательные постоянные gk в этой системе являются известными, а µ - коэффициент трения. Значения параметров bi , aij , gk в (24) с учетом (22), (23) определяются следующими соотношениями: 1 1 ˜b1 (s)ds, b ˜b2 (s)ds, 3k-1 = mes(ωk ) ωk mes(ωk ) ωk 1 1 ˜b3 (s)ds, b3k = gk = g(s)ds, mes(ωk ) ωk mes(ωk ) ωk 1 K3-d 3-e (t, s)ds dt + εδ3k-d 3j-e , a3k-d 3j-e = mes(ωk ) ωk ωj k, j = 1, n2 , d, e = 0, 2. b3k-2 = (25) Здесь b = (˜b1 , ˜b2 , ˜b3 ), δij - символ Кронекера, а положительное число ε - параметр регуляризации. Предложенный способ построения аппроксимирующих последовательностей {An }, {bn }, {gn } гарантирует, что решение pn ∈ L32 (Ω) операторного уравнения (19) (существование и единственность которого в L32 (Ω) при достаточно большом n вытекает из теоремы 4) непременно должно содержаться в конечномерном подпространстве Hn пространства L32 (Ω). Это означает, что при выполнении всех условий теоремы 4 (а также при полной непрерывности линейных интегральных операторов Aij : L2 (Ω) → L2 (Ω)) система уравнений (24) является совместной, если только значение n достаточно велико. Таким образом, искомое решение pn уравнения (19) при достаточно большом n может быть восстановлено по найденным значениям неизвестных x1 , x2 , . . . , x3n2 , удовлетворяющих системе уравнений (24). Для многих контактных задач теории упругости используется аппроксимация взаимодействующих тел упругими полупространствами [10]. В этом случае для последовательности {An }, определяемой формулами (22), (23), можно вместо формул (23) использовать следующие соотношения [9]:  n n  Kd e (si , sj ), если t ∈ ωi , s ∈ ωj , i = j; (n) (26) Kd e (t, s) = 1  Kd e (sni , s)ds, если t, s ∈ ωi ,  mes(ωi ) ωi где i, j = 1, n2 , d, e = 1, 3. В соотношениях (26) sni - центр квадрата ωi , а ядра Kij (t, s) интегральных операторов Aij определяются в соответствии с формулами Буссинеска-Черрути [10]. Таким образом, в случае аппроксимации взаимодействующих тел упругими полупространствами коэффициенты aij в системе (24) можно вычислить без использования соотношений (25) по следующим простым формулам: a3k-2 3k-1 = a3k-2 3k = a3k-1 3k-2 = a3k-1 3k = a3k 3k-2 = a3k 3k-1 = 0, a3k-d 3j-e = mes(ωj )K3-d 3-e (snk , snj ), K3-j 3-j (snk , s)ds, a3k-d 3k-d = ε + ωk ∀k, j = 1, n2 (k = j), 320 ∀d, e = 0, 2. (27) Метод нелинейных граничных интегральных уравнений. . . Отметим, что при аппроксимации взаимодействующих тел упругими полупространствами линейные интегральные операторы Aij : L2 (Ω) → L2 (Ω) являются вполне непрерывными [15], а порождаемый ими вполне непрерывный оператор A : L32 (Ω) → L32 (Ω) вида (8) удовлетворяет всем условиям теорем 1-4. 6. Итерационный процесс для решения дискретного аналога регуляризированного уравнения. Из всех доказанных выше утверждений следует, что нахождение приближенного решения операторного уравнения (11), соответствующего определенному шагу нагружения рассматриваемой квазистатической контактной задачи, сводится к решению системы нелинейных скалярных уравнений (24) при достаточно малом ε > 0 и достаточно большом натуральном n. Для приближенного решения системы (24) можно предложить следующий итерационный процесс: (0) (0) (0) 2 (x1 , x2 , . . . , x3n2 ) ∈ R3n ; 3n2 (m+1) x3k-2 =h (m) x3k-2 -E (m) - b3k-2 (m) - b3k-1 , a3k-2 j xj ; j=1 3n2 (m+1) (m) x3k-1 = q x3k-1 - E a3k-1 j xj j=1 3n2 (m) x3k (m) -E a3k j xj - b3k , µgk ; (28) j=1 3n2 (m+1) x3k (m) (m) = q x3k - E a3k j xj - b3k , j=1 3n2 (m) x3k-1 (m) -E a3k-1 j xj - b3k-1 , µgk ; j=1 k = 1, n2 , m = 0, 1, 2, . . . Сходимость итерационного процесса (28) к решению системы (24) будет иметь (0) (0) (0) место независимо от выбора начального приближения (x1 , x2 , . . . , x3n2 ) 2 в 3n2 -мерном евклидовом пространстве R3n при условии, что матрица A0 = 2 2 ,3n = (aij )3n i=1, j=1 является симметрической и положительно определенной, а входящая в (28) положительная константа E удовлетворяет неравенству 3n2 E< 1 i 3n2 -1 |aij | max . j=1 Этот факт следует из принципа сжимающих отображений. 7. Численные результаты. Предложенным методом было получено численное решение пространственной квазистатической контактной задачи о взаимодействии упругого шара радиуса 0.3 м с упругим полупространством. Процесс нагружения тел осуществлялся в соответствии со следующим законом 321 С т р е л я е в Ю. М. изменения жёстких смещений ∆1i , ∆2i , ∆3i шара, выраженных в метрах:   -0.000002 · i, если 1 i 50; ∆1i = (29)  -0.0002 + 0.000002 · i, если 51 i 65; ∆2i = ∆3i = 0, если 1 i 65. График зависимости абсолютной величины нормального жесткого смещения шара |∆1i | от номера шага нагружения i представлен на рис. 2. В соотношениях (29) величины ∆2i и ∆3i есть касательные составляющие жесткого смещения шара в направлении осей Ox и Oy соответственно. Рис. 2. Кривая нагружения [Figure 2. A loading curve] В расчетах использовались следующие модули Юнга E1 , E2 и коэффициенты Пуассона ν1 , ν2 шара и полупространства: E1 = 3 · 105 МПа, E2 = = 105 МПа, ν1 = ν2 = 0.2. Коэффициент трения µ = 0.12375. Значение безразмерного параметра β, определяемого согласно [10] формулой β= 1 E1 (1 - 2ν2 ) (1 + ν2 ) - E2 (1 - 2ν1 ) (1 + ν1 ) , 2 E1 1 - ν22 + E2 1 - ν12 составляет для рассматриваемой контактной задачи 0.1875. Отношение µ/β = = 0.66. При дискретизации зона контакта разбивалась на 41 × 41 = 1681 квадратных граничных элемента одинаковой площади (n = 41, mes(ωi ) = = 7.84 · 10-8 м2 ). Для вычисления элементов akj матрицы податливости использовались соотношения (27), в которых функции Kjm (s, s ) задавались при помощи формул Буссинеска-Черрути, а значение ε определялось соотношением ε = 10-3 · K11 (snj , s)ds. ωk С помощью итерационного процесса (28) в настоящей работе были получены два численных решения контактной задачи, которые соответствуют различным этапам процесса нагружения. Первое решение было получено на 322 Метод нелинейных граничных интегральных уравнений. . . тридцать пятом шаге нагружения, что соответствует монотонному нормальному внедрению шара в полупространство на 0.00007 м (отрезок OA, рис. 2). Второе решение, полученное на последнем, шестьдесят пятом, шаге нагружения (ломаная OBC, рис. 2), соответствует первоначальному вдавливанию шара в полупространство на 0.0001 м (отрезок OB, рис. 2) и последующей разгрузке (отрезок BC, рис. 2). Отметим, что одноименные значения относительных жёстких смещений тел на тридцать пятом шаге совпали с их значениями на шестьдесят пятом шаге процесса нагружения тел. Полученные из первого решения значения касательных напряжений τzx , действующих в точках поверхности шара, соприкасающихся с полупространством по оси Ox, сопоставлялись с результатами решения соответствующей осесимметричной контактной задачи, приведенными в работе [10] для случая µ/β = 0.66. В таблице показана зависимость безразмерной величины τzx /(βpmax ) от безразмерного параметра x/a. Первый столбец таблицы соответствует решению [10], второй - численному решению, полученному на тридцать пятом шаге процесса нагружения (29), третий - численному решению, полученному на шестьдесят пятом шаге этого процесса. Здесь символ pmax обозначает максимальное значение контактного давления, а символ a - радиус площадки контакта. Данные, приведенные в первом и втором столбцах таблицы, показывают хорошее соответствие решения [10] и численного решения, полученного при монотонном нагружении. Максимальная относительная погрешность в этом случае составляет приблизительно 6 % по отношению к наибольшему значению величины τzx /(βpmax ), найденному согласно [10], что свидетельствует о корректности предложенного метода решения контактной задачи. Результаты, приведенные в третьем столбце таблицы, означают, что при разгрузке наблюдается уменьшение достигнутых значений касательных напряжений и смена знака этих значений вблизи границы пятна контакта. Зависимость безразмерных касательных напряжений от безразмерного расстояния [The dependence of the dimensionless shear stresses on the dimensionless distance] Solution by the Eqs. (28), (29) Solution by the Eqs. (28), (29) x/a Solution from [10] at the step 35 at the step 65 0.00 0.06 0.12 0.18 0.24 0.30 0.36 0.42 0.48 0.54 0.60 0.66 0.72 0.78 0.84 0.90 0.96 0.00 0.10 0.18 0.25 0.30 0.35 0.39 0.42 0.45 0.47 0.48 0.48 0.45 0.40 0.33 0.23 0.03 0.00 0.07 0.17 0.22 0.29 0.35 0.39 0.42 0.45 0.47 0.47 0.46 0.42 0.37 0.30 0.22 0.05 0.00 0.07 0.17 0.23 0.30 0.35 0.40 0.43 0.46 0.48 0.50 0.48 0.44 0.40 0.32 0.12 -0.14 323 С т р е л я е в Ю. М. На рис. 3 показаны конфигурации пятен контакта и зон сцепления, соответствующие полученным численным решениям. Здесь внешний замкнутый контур есть граница пятна контакта, заштрихованная область - зона сцепления, незаштрихованная область - зона проскальзывания. Рис. 3, a соответствует решению, полученному на тридцать пятом шаге, рис. 3, b - решению, полученному на шестьдесят пятом шаге процесса нагружения (29). a b Рис. 3. Границы областей контакта и зон сцепления шара и полупространства при монотонном (a) и немонотонном (b) нагружении [Figure 3. The boundaries of the contact areas and the adhered regions of sphere and half-space under monotonic (a) and nonmonotonic (b) loading] Анализ представленных на рис. 3 результатов показал, что качественные характеристики контакта тел невозможно установить по конечным значениям жёстких смещений этих тел, не учитывая, каким образом эти конечные значения были достигнуты в процессе нагружения. Полученные данные также показывают, что после прохождения этапа увеличения нагрузки до максимального значения (отрезок AB, рис. 2) и последующей разгрузки (отрезок BC, рис. 2) граница центрального участка зоны сцепления несколько сместилась к центру и одновременно вблизи границы пятна контакта образовалась вторая кольцевая зона сцепления. Отметим, что подобную конфигурацию зон сцепления на контактном пятне получил J. R. Turner [16] при исследовании случая разгрузки в задаче о взаимодействии цилиндрического штампа, имеющего плоскую подошву, с упругим полупространством. Выводы. Разработан метод получения приближённого решения трёхмерной квазистатической контактной задачи о взаимодействии двух линейноупругих тел с учетом кулонова трения между ними при неизвестной заранее площадке контакта и неизвестной границе раздела зон проскальзывания и сцепления на этой площадке. Анализ полученных численных результатов для задачи о вдавливании упругого шара в упругое полупространство свидетельствует о корректности предложенного метода и его применимости к решению контактных задач рассматриваемого класса. 324 Метод нелинейных граничных интегральных уравнений. . .

About the authors

Yurii M Streliaiev

Zaporizhzhya National University

Email: strelkiny@gmail.com
(strelkiny@gmail.com), Senior Lecturer, Dept. of Mathematical Analysis 66, Zhukovskogo st., Zaporizhzhya, 69600, Ukraine

References

Supplementary files