Известия РАН. Механика твердого тела, 2023, № 2, стр. 70-89

1. Введение. При контактном взаимодействии твердых тел область фактического контакта, как правило, дискретна. Размеры и положение пятен фактического контакта зависят от условий контактного взаимодействия, механических характеристик тел и их поверхностной микроструктуры. Для описания дискретного (множественного) контакта твердых тел предложены различные математические модели, учитывающие параметры макро- и микрогеометрии реальных поверхностей и условия их контактного взаимодействия [1–4]. Большинство этих моделей основано на классической теории контактного взаимодействия деформируемых тел [5, 6]. Подробный обзор современного состояния исследований в области механики дискретного контакта, включая основные подходы к постановке задач, методы аналитического и численного решения, конкретные результаты и области их практического использования, приведен в статье [7].

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

Задачи дискретного контакта с односторонними связями являются нелинейными вследствие наличия в постановке граничных условий в виде неравенств. Для их численного решения требуется применение итерационных алгоритмов. Наиболее распространенный подход к построению вычислительных алгоритмов решения контактных задач с односторонними связями состоит в использовании вариационных методов [8–10]. Вариационные формулировки также применяются для исследования проблемы существования, единственности и регулярности решения контактных задач [11–13].

Используя математический аппарат преобразования вариационных задач [8], можно получить семейство вариационных формулировок. С вычислительной точки зрения наиболее эффективным является применение граничных вариационных формулировок с использованием оператора Пуанкаре–Стеклова, отображающего на части границы упругого тела, по которой возможен контакт, поверхностные силы в перемещения поверхности. В этом случае для дискретизации контактной задачи применяется метод граничных элементов, причем дискретизации подлежит не вся граница упругого тела, а лишь ее часть – область возможного контакта.

Для построения оператора Пуанкаре–Стеклова требуется решить линейную краевую задачу с заданными силами на поверхности возможного контакта. Как известно [14], при решении краевых задач для упругого слоя может быть эффективно использовано двумерное интегральное преобразование Фурье. В ряде случаев удается получить аналитические выражения для трансформант Фурье искомых функций, однако выполнить аналитически обратное преобразование Фурье этих двумерных трансформант весьма затруднительно, поэтому используют приближенные методы. Асимптотические методы представлены в монографиях [15, 16]. Численные методы, как правило, основаны на двумерном дискретном преобразовании Фурье (ДПФ), при реализации которого применяются алгоритмы быстрого преобразования Фурье (БПФ) [17, 18]. Подробный обзор современного состояния исследований и анализ различных аспектов применения ДПФ и БПФ в контактных задачах приведен в [19, 20].

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

Впервые алгоритм решения задач квадратичного программирования на основе МСГ был предложен в работе [21]. Идея алгоритма состояла в использовании стратегии рабочего списка активных ограничений-неравенств (активного набора) и применении МСГ для нахождения минимума на текущем рабочем подпространстве. Была доказана теоретическая сходимость алгоритма за конечное число шагов. При практическом использовании алгоритма оказалось, что скорость сходимости итерационного процесса существенно зависит от количества и вида ограничений. Поэтому впоследствии были предложены различные модификации базового алгоритма, адаптированные для решения конкретных классов прикладных задач. Описание и анализ алгоритмов на основе МСГ для решения контактных задач с односторонними связями можно найти, например, в работах [22, 23].

Один из наиболее известных алгоритмов, приведенный в часто цитируемой статье [24], отличается от базового алгоритма [21] учетом специфики множества ограничений решаемой задачи квадратичного программирования, что позволило значительно повысить скорость работы алгоритма. Однако алгоритм [24] имеет существенное ограничение. В нем учитывается только одно уравнение равновесия штампа (для сил) и, следовательно, только поступательное движение штампа. Это не позволяет использовать алгоритм для решения контактных задач при внецентренном нагружении штампа, когда необходимо учитывать его повороты.

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

В настоящей работе алгоритм [25] обобщен для пространственных контактных задач. Предложено линейное преобразование переменных, позволяющее учесть все три уравнения равновесия штампа и рассматривать случаи внецентренного нагружения штампа в пространственных контактных задачах для упругого слоя. Отметим также, что в отличие от работы [25] при аппроксимации оператора Пуанкаре–Стеклова используется ДПФ, для численной реализации которого применяются алгоритмы БПФ. Это позволило существенно сократить вычислительные затраты при решении пространственных контактных задач.

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

2. Постановка задачи. Пусть в неподвижной прямоугольной системе координат $O{{x}_{1}}{{x}_{2}}{{x}_{3}}$ невесомый однородный изотропный упругий слой занимает область Ω = = $\{ x = ({{x}_{1}},{{x}_{2}},{{x}_{3}}) \in {{\mathbb{R}}^{3}}:\;\left| {{{x}_{1}}} \right| \leqslant \infty ,\;\left| {{{x}_{2}}} \right| \leqslant \infty ,\;0 \leqslant {{x}_{3}} \leqslant h\} $. Границу слоя ${{x}_{3}} \equiv 0$ будем обозначать ${{\Gamma }_{0}}$, а границу ${{x}_{3}} \equiv h$ – через ${{\Gamma }_{1}}$. Далее под $u(x),\;\varepsilon (x),\;\sigma (x)$ будем понимать соответственно вектор перемещений и тензоры деформаций и напряжений в точке $x \in \Omega $. Предполагается, что деформации малы, а напряжения в недеформированном состоянии отсутствуют. Напряженно-деформированное состояние слоя описывается системой уравнений [14, 15]:

где ${\text{def}} \equiv 1{\text{/}}2({\text{grad}} + {\text{gra}}{{{\text{d}}}^{T}})$, $S$ – тензор модулей упругости.

По границе ${{\Gamma }_{0}}$ слой соединен с недеформируемым основанием. В случае идеального контакта граничные условия имеют вид

В слой вдавливается гладкий жесткий штамп, основание которого имеет поверхностный рельеф. Часть границы ${{\Gamma }_{1}}$, по которой возможен контакт слоя со штампом, обозначается Γ_p. Положение и предельные размеры Γ_p, т.е. номинальная область контакта, задаются априори, исходя из геометрических соображений. Без потери общности будем полагать, что область Γ_p представляет собой прямоугольник размеров ${{b}_{1}} \times {{b}_{2}}$ со сторонами, параллельными осям координат $O{{x}_{1}}$ и Ox₂ соответственно. Отметим, что при вдавливании штампа с поверхностным рельефом номинальная область контакта Γ_p включает в себя множество отдельных пятен фактического контакта, положение и форма которых заранее неизвестны и подлежат определению в процессе решения задачи.

Макроформа основания штампа и его поверхностный рельеф описываются функцией $\Phi {\kern 1pt} ({{x}_{1}},{{x}_{2}})$, значение которой в точке $x \in {{\Gamma }_{p}}$ равно расстоянию от этой точки до поверхности штампа, измеренному вдоль направления внешней нормали к границе Γ_p. Расстояние $\Phi {\kern 1pt} ({{x}_{1}},{{x}_{2}})$ отсчитывается по отношению к недеформированному состоянию слоя. Для определенности будем полагать $\mathop {min}\limits_{x \in {{\Gamma }_{p}}} \Phi {\kern 1pt} ({{x}_{1}},{{x}_{2}}) = 0$. Для штампа с поверхностным рельефом функция $\Phi {\kern 1pt} ({{x}_{1}},{{x}_{2}})$ является мультимодальной (многоэкстремальной). Положение штампа определяется векторами перемещений $\delta = ({{\delta }_{1}},{{\delta }_{2}},{{\delta }_{3}})$ и углов поворота $\varphi = ({{\varphi }_{1}},{{\varphi }_{2}},{{\varphi }_{3}})$. Главный вектор $F = ({{F}_{1}},{{F}_{2}},{{F}_{3}})$ и главный момент M = $({{M}_{1}},{{M}_{2}},{{M}_{3}})$ внешних сил, приложенных к штампу, считаются заданными. В качестве центра приведения выбирается точка ${{x}^{c}} = (x_{1}^{c},x_{2}^{c},x_{3}^{c})$. Далее рассматривается задача нормального контакта слоя со штампом, поэтому будем полагать ${{F}_{1}} = {{F}_{2}} = 0$, ${{M}_{3}} = 0$, ${{\delta }_{1}} = {{\delta }_{2}} = 0$ и ${{\varphi }_{3}} = 0$.

Контактное взаимодействие упругого слоя с жестким штампом описывается условиями одностороннего гладкого контакта:

Остальная часть границы ${{\Gamma }_{1}}$ свободна от внешних нагрузок:

Уравнения равновесия жесткого штампа имеют вид:

Отметим, что соотношения (2.5), по существу, представляют собой нелокальные граничные условия.

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

– компонента главного вектора ${{F}_{3}}$ ограничена по абсолютной величине и ${{F}_{3}} < 0$;

– значения силы ${{F}_{3}}$ и моментов ${{M}_{1}}$ и ${{M}_{2}}$ согласованы между собой таким образом, что существует распределение нормальных напряжений ${{\sigma }_{{33}}} \leqslant 0$ на Γ_p, удовлетворяющее уравнениям равновесия штампа (2.5).

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

которое вполне замыкает постановку задачи и определяет поведение решения на бесконечности [15].

Задача (в дифференциальной постановке) состоит в определении полей перемещений $u(x)$, деформаций $\varepsilon (x)$ и напряжений $\sigma (x)$, удовлетворяющих уравнениям (2.1), граничным условиям (2.2)–(2.4), уравнениям равновесия штампа (2.5) и условию (2.6). Также необходимо найти смещение ${{\delta }_{3}}$ и повороты ${{\varphi }_{1}}$ и ${{\varphi }_{2}}$ жесткого штампа. Подчеркнем, что в рассматриваемой задаче дискретного контакта априори задается лишь предельно допустимая (номинальная) область контакта Γ_p, положение и форма пятен фактического контакта заранее неизвестны и подлежат определению.

Отметим, что более общие постановки задач дискретного контакта для упругого слоя с заданными массовыми силами и неоднородными смешанными граничными условиями на ${{\Gamma }_{1}}{{\backslash }}{{\Gamma }_{p}}$ можно рассматривать как контактные задачи при наличии пригрузки вне штампа [1, 5]. Следуя общему подходу, несложно показать, что эти задачи могут быть сведены к рассматриваемой путем модификации функции, описывающей форму основания штампа

где $\tilde {u} = ({{\tilde {u}}_{1}},{{\tilde {u}}_{2}},{{\tilde {u}}_{3}})$ – решение краевой задачи для упругого слоя с заданными массовыми силами и неоднородными смешанными граничными условиями на ${{\Gamma }_{1}}{{\backslash }}{{\Gamma }_{p}}$ при отсутствии штампа, т.е. когда

3. Функциональные пространства. Решение $u(x)$ задачи (2.1)–(2.6), принадлежащее классу функций ${{[{{C}^{2}}(\Omega )]}^{3}} \cap {{[{{C}^{1}}(\bar {\Omega })]}^{3}}$, является классическим решением. Переход к вариационной формулировке задачи позволяет определить обобщенное решение. Введем необходимые для этого функциональные пространства [15, 26, 27]. Обозначим

где ${{H}^{1}}(\Omega ) = {{[{{H}^{1}}(\Omega )]}^{3}},$ ${{L}_{2}}(\Omega ) = {{\left[ {{{L}_{2}}(\Omega )} \right]}^{3}}$ – прямые произведения пространств Соболева. Пространство ${{H}^{1}}(div;\Omega )$ является гильбертовым со скалярным произведением

Для однородного изотропного упругого слоя несложно получить неравенство

с положительной постоянной c₁, не зависящей от $v$. Следовательно, элементы пространства ${{H}^{1}}(div;\Omega )$ удовлетворяют условию конечности потенциальной энергии деформации (2.6).

Выделим в ${{H}^{1}}(div;\Omega )$ подпространство

Нетрудно видеть, что

Из результатов [11] с учетом (3.5) следует, что существует постоянная ${{c}_{2}} > 0$, зависящая лишь от $\Omega $, такая, что

Неравенство (3.6) представляет собой первое неравенство Корна.

Используя далее результаты [15, 28], можно показать, что существуют линейные непрерывные операторы следа

определяющие для полей перемещений $v \in V$ соответственно векторные функции перемещений ${{\gamma }_{u}}v = ({{{v}}_{1}},{{{v}}_{2}},{{{v}}_{3}})$ и напряжений ${{\gamma }_{t}}v = ({{\sigma }_{{31}}},{{\sigma }_{{32}}},{{\sigma }_{{33}}})$ на ${{\Gamma }_{1}}$ такие, что где c₃, c₄ – положительные постоянные, не зависящие от v.

Обратно, существует линейный непрерывный оператор

такой, что для заданного $t \in {{[{{H}^{{ - 1{\text{/}}2}}}({{\Gamma }_{1}})]}^{3}}$ можно построить функцию $v \in V$, обладающую свойством ${{\gamma }_{t}}{\kern 1pt} v = t$.

Введем пространство сужений на Γ_p функций из ${{H}^{{1{\text{/}}2}}}({{\Gamma }_{1}})$

и оснастим его нормой

Имеет место вложение пространств ${{H}^{{1{\text{/}}2}}}({{\Gamma }_{p}}) \subset {{L}_{2}}({{\Gamma }_{p}})$ [29]. Двойственным к пространству ${{H}^{{1{\text{/}}2}}}({{\Gamma }_{p}})$ является пространство ${{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$ – пополнение ${{L}_{2}}({{\Gamma }_{p}})$ по норме

т.е. отношение двойственности ${}_{{{{{\tilde {H}}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})}}{{\langle \cdot , \cdot \rangle }_{{{{H}^{{1{\text{/}}2}}}({{\Gamma }_{p}})}}}$ порождается продолжением стандартного скалярного произведения в ${{L}_{2}}({{\Gamma }_{p}})$ и имеет место плотное вложение пространств [29]

Для упрощения обозначений указанное отношение двойственности далее будем обозначать как $\langle \cdot , \cdot \rangle $.

Заметим, что тривиальное продолжение функции $p \in {{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$ на Γ₁ принадлежит пространству ${{H}^{{ - 1{\text{/}}2}}}({{\Gamma }_{1}})$, т.е. существует функция $p{\kern 1pt} * \in {{H}^{{ - 1{\text{/}}2}}}({{\Gamma }_{1}})$ такая, что supp$p{\kern 1pt} * \subset {{\bar {\Gamma }}_{p}}$, p является сужением p* на Γ_p и выполняется равенство [28]

Следовательно, с учетом (3.9) существует линейный непрерывный оператор

такой, что для заданного $p \in {{[{{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})]}^{3}}$ можно построить функцию $v \in V$, обладающую свойством ${{\left. {{{\gamma }_{t}}{\kern 1pt} v} \right|}_{{{{\Gamma }_{p}}}}} = p$.

4. Интегральное представление решения. Рассмотрим вспомогательную краевую задачу для упругого слоя: найти поле перемещений u, удовлетворяющее уравнениям (2.1), граничным условиям (2.2) и (2.4), а также условиям

В [15] показано, что для любого $q \in {{L}_{2}}({{\Gamma }_{p}})$ данная задача имеет единственное решение $u$ и выполняется интегральное тождество (обобщенная формула Грина)

Из (4.2) с учетом (3.6) следует оценка

с положительной постоянной c₅, не зависящей от q. Следовательно, с помощью решения вспомогательной краевой задачи можно определить линейный непрерывный оператор ${{G}_{p}}:q \mapsto u$ из ${{L}_{2}}({{\Gamma }_{p}})$ в V.

Решение вспомогательной краевой задачи для $q \in {{L}_{2}}({{\Gamma }_{p}})$ имеет вид [16]:

Выражение для ядра $U(\alpha ,\beta ,{{x}_{3}})$ приведено, например, в [16]. Учитывая плотность вложения (3.13), оператор G_p можно продолжить по непрерывности на ${{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$ и с учетом (3.15) рассматривать как оператор, действующий из ${{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$ в $V$.

Таким образом, использование интегрального представления (4.4)–(4.5) в предположении, что $\Phi \in {{H}^{{1{\text{/}}2}}}({{\Gamma }_{p}})$, позволяет свести решение контактной задачи (2.1)–(2.6) к нахождению нормальных напряжений ${{\sigma }_{{33}}} = q \in {{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$, удовлетворяющих следующей системе уравнений и неравенств:

Перемещения $u = {{G}_{p}}q$, соответствующие решению $q \in {{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$ системы (4.6)–(4.9), будем называть обобщенным решением контактной задачи (2.1)–(2.6).

5. Оператор Пуанкаре–Стеклова. Введем оператор Пуанкаре–Стеклова ${{G}_{s}}:q \mapsto w$, отображающий посредством решения (4.4)–(4.5) нормальные напряжения q(x₁, x₂) ≡ ≡ ${{\sigma }_{{33}}}({{x}_{1}},{{x}_{2}},h)$ на части Γ_p границы упругого слоя в нормальные перемещения w(x₁, x₂) ≡ ≡ ${{u}_{2}}({{x}_{1}},{{x}_{2}},h)$ на Γ_p.

Соответствующее решению (4.4)–(4.5) выражение для оператора Пуанкаре–Стеклова ${{G}_{s}}:q \mapsto w$ для $p \in {{L}_{2}}({{\Gamma }_{p}})$ имеет вид [16]:

Отметим, что ядро ${{K}_{s}}( \cdot , \cdot )$ интегрального представления (5.1) является разностным по обоим аргументам [16]. Учитывая плотность вложения (3.13), оператор ${{G}_{s}}$ можно продолжить по непрерывности на ${{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$ и с учетом (3.7), (3.10), (3.15) рассматривать как оператор, действующий из ${{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$ в ${{H}^{{1{\text{/}}2}}}({{\Gamma }_{p}})$.

Введем на ${{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}}) \times {{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$ непрерывную билинейную форму

Учитывая плотность вложения (3.17), из (4.2) получим, что

где $u = {{G}_{p}}q$, $v = {{G}_{p}}p$.

Для изотропного упругого слоя билинейная форма $g( \cdot , \cdot )$ является симметричной. Используя (3.6), далее получим неравенство

из которого с учетом (3.8) и (3.14) следует оценка с положительной постоянной c₆, не зависящей от q. В результате, билинейная форма $g( \cdot , \cdot )$ и оператор Пуанкаре–Стеклова ${{G}_{s}}$ являются положительно определенными на ${{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$.

6. Вариационная формулировка задачи. Для решения системы (4.6)–(4.9), содержащей уравнения и неравенства, используется вариационный подход [8–13]. Образуем множество статически допустимых нормальных напряжений на ${{\Gamma }_{p}}$

Нетрудно видеть, что множество Σ является замкнутым выпуклым множеством в ${{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$. Кроме того, в соответствии со сделанными при постановке задачи предположениями это множество является непустым.

Пусть $q \in \Sigma $ – решение системы (4.6)–(4.9). Тогда для произвольного $p \in \Sigma $ получим следующую оценку

Из (6.2) следует, что искомые нормальные напряжения $q \in \Sigma $ удовлетворяют граничному вариационному неравенству

где $\phi ( \cdot )$ – определенная на ${{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$ линейная форма

Отметим, что неравенство (6.3) не содержит неизвестных смещения ${{\delta }_{3}}$ и поворотов ${{\varphi }_{1}}$ и ${{\varphi }_{2}}$ жесткого штампа. Как нетрудно видеть из преобразований (6.2), это является следствием того, что элементы множества статически допустимых нормальных напряжений Σ удовлетворяют уравнениям равновесия штампа (4.9).

Используя известные приемы [30, 31], можно доказать обратное утверждение: решение q вариационного неравенства (6.3) удовлетворяет системе уравнений и неравенств (4.6)–(4.9).

Учитывая, что билинейная форма $g(p,q)$ является положительно определенной, а множество Σ – непустым замкнутым выпуклым множеством в ${{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$, вариационное неравенство (6.3) эквивалентно задаче минимизации граничного функционала: найти $q \in \Sigma $ такой, что

Кроме того, решение вариационного неравенства (6.3) и задачи минимизации (6.5) существует и единственно [8, 12, 13].

Отыскав нормальные напряжения $q \equiv {{\sigma }_{{33}}}$ на Γ_p как решение задачи (6.5), можно определить напряженно-деформированное состояние всего слоя, используя (4.4)–(4.5).

7. Аппроксимация задачи. Выбор метода аппроксимации задачи минимизации (6.5) определяется, в первую очередь, возможностью построения эффективной с вычислительной точки зрения аппроксимации оператора Пуанкаре–Стеклова ${{G}_{s}}$, ядро интегрального представления (5.1) которого имеет вид [16]

где

Аналитическое выражение для интеграла (7.1) неизвестно. Применение кубатурных формул для вычисления этого несобственного двойного интеграла представляет собой трудоемкую вычислительную задачу. Альтернативный подход состоит в следующем. Нетрудно видеть, что ${{\tilde {K}}_{s}}(\alpha ,\beta )$ является двумерной косинус-трансформантой Фурье функции ${{K}_{s}}({{z}_{1}},{{z}_{2}})$. Подвергнем (5.1) двумерному преобразованию Фурье и в соответствии с теоремой о свертке [32] получим

где $\tilde {w}(\alpha ,\beta ),$ $\tilde {q}(\alpha ,\beta )$ – двумерные трансформанты Фурье соответственно нормальных перемещений $w({{x}_{1}},{{x}_{2}})$ на ${{\Gamma }_{1}}$ и тривиального продолжения на ${{\Gamma }_{1}}$ нормальных напряжений $q({{x}_{1}},{{x}_{2}})$. Следовательно, вычисление оператора Пуанкаре–Стеклова ${{G}_{s}}$ сводится к выполнению пары (прямого и обратного) двумерных преобразований Фурье и перемножению спектров (7.3).

В численном анализе использование преобразования Фурье наиболее эффективно в случае периодических функций благодаря дискретности их спектра и, как следствие, переходе от непрерывного преобразования к дискретному [17, 18]. Поэтому основная идея используемого ниже подхода состоит в аппроксимации искомых нормальных напряжений на ${{\Gamma }_{p}}$ сеточными функциями, периодическими на всей границе ${{\Gamma }_{1}}$, и применении алгоритмов БПФ. Для уменьшения возникающей при этом так называемой ошибки периодичности вводится расширенная вычислительная область ${{\Gamma }_{c}}$ [19, 20].

Выберем прямоугольную вычислительную область ${{\Gamma }_{c}} \subset {{\Gamma }_{1}}$ размеров ${{a}_{1}} \times {{a}_{2}}$ со сторонами, параллельными осям координат $O{{x}_{1}}$ и $O{{x}_{2}}$, так, чтобы выполнялись условия

где ${{\chi }_{i}}$ – коэффициенты расширения вычислительной области. Далее аналогично (3.10)–(3.12) введем пространство ${{H}^{{1{\text{/}}2}}}({{\Gamma }_{c}})$ и двойственное ему пространство ${{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{c}})$. Используя результаты [29], можно показать, что сужения функций из ${{H}^{{1{\text{/}}2}}}({{\Gamma }_{c}})$ на ${{\Gamma }_{p}}$ принадлежат пространству ${{H}^{{1{\text{/}}2}}}({{\Gamma }_{p}})$, а тривиальное продолжение ${{q}_{0}}$ функции $q \in {{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$ на Γ_c принадлежит пространству ${{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{c}})$. Обозначим через ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{H} }^{{ - 1{\text{/}}2}}}({{\Gamma }_{c}})$ подпространство функций из ${{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{c}})$, имеющих носитель на ${{\bar {\Gamma }}_{p}}$. Учитывая (7.4), аналогично [33] продолжим периодически функции из ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{H} }^{{ - 1{\text{/}}2}}}({{\Gamma }_{c}})$ на всю границу Γ₁ с периодом Γ_c и образуем пространство $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{H} _{{per}}^{{ - 1{\text{/}}2}}({{\Gamma }_{1}})$ Γ_c-периодических функций. Далее функции из ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{H} }^{{ - 1{\text{/}}2}}}({{\Gamma }_{c}})$ будем отождествлять с их периодическими продолжениями на Γ₁. Тем самым определим оператор $\Pi :{{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}}) \to \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{H} _{{per}}^{{ - 1{\text{/}}2}}({{\Gamma }_{1}})$, ставящий в соответствие функции из $q \in {{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{p}})$ периодическое продолжение ${{q}_{{per}}} \in \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{H} _{{per}}^{{ - 1{\text{/}}2}}({{\Gamma }_{1}})$ на Γ₁ ее тривиального продолжения ${{q}_{0}} \in {{\tilde {H}}^{{ - 1{\text{/}}2}}}({{\Gamma }_{c}})$ на Γ_c.

Обозначим через ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{L} }_{2}}({{\Gamma }_{c}})$ подпространство функций из ${{L}_{2}}({{\Gamma }_{c}})$, имеющих носитель на ${{\bar {\Gamma }}_{p}}$. Используя (3.13), можно показать, что имеет место плотное вложение пространств ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{L} }_{2}}({{\Gamma }_{c}}) \subset {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{H} }^{{ - 1{\text{/}}2}}}({{\Gamma }_{c}})$. Функции из ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{L} }_{2}}({{\Gamma }_{c}})$ можно рассматривать как Γ_c-периодические, следовательно система функций $exp\left[ {2\pi i\left( {{{l}_{1}}{{x}_{1}}{\text{/}}{{a}_{1}} + {{l}_{2}}{{x}_{2}}{\text{/}}{{a}_{2}}} \right)} \right],$ ${{l}_{1}},{{l}_{2}} = 0, \pm 1, \pm 2, \ldots ,$ является плотной в ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{L} }_{2}}({{\Gamma }_{c}})$ [34]. Поэтому $q \in {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{H} }^{{ - 1{\text{/}}2}}}({{\Gamma }_{c}})$ можно аппроксимировать с заданной точностью конечным отрезком ряда Фурье

Для вычисления коэффициентов Фурье $c{\kern 1pt} {\kern 1pt} [{{l}_{1}},{{l}_{2}}] \in \mathbb{Z}$ с помощью ДПФ введем сеточные функции. Построим на ${{\Gamma }_{c}}$ регулярную (равномерную) сетку ${{\Upsilon }_{c}}$, состоящую из ${{N}_{c}} = {{N}_{{c1}}} \cdot {{N}_{{c2}}}$ одноузловых прямоугольных граничных элементов нулевого порядка размеров ${{d}_{1}} \times {{d}_{2}}$, где ${{d}_{i}} = {{a}_{i}}{\text{/}}{{N}_{{ci}}},\;i = 1,2$. Узлы ${{\Upsilon }_{c}}$ обозначим ${{s}_{{nm}}},\;n = \overline {1,{{N}_{{c1}}}} ,\;m = \overline {1,{{N}_{{c2}}}} $. Часть ${{\Upsilon }_{c}}$, расположенную на Γ_p и содержащую $N = {{N}_{{p1}}} \cdot {{N}_{{p2}}}$ элементов, обозначим ${{\Upsilon }_{p}}$. Из (7.4) следует, что ${{N}_{{pi}}} = {{N}_{{ci}}}{\text{/}}{{\chi }_{i}},\;i = 1,2$.

Далее векторы и матрицы будем обозначать большими латинскими буквами, а их элементы – соответствующими малыми латинскими буквами.

Введем на ${{\Upsilon }_{c}}$ сеточные функции (матрицы) нормальных напряжений и перемещений

и соответствующие им векторы где $vec( \cdot )$ – оператор векторизации матрицы, преобразующий ее в вектор путем расположения ее столбцов друг под другом.

Введем также редукции сеточных функций на ${{\Upsilon }_{p}}$

где ${{R}_{{pc}}}$ – прямоугольная матрица размеров $N \times {{N}_{c}}$, в каждой строке которой имеется ровно один ненулевой элемент, равный 1.

Обратные операции тривиального продолжения (дополнения нулевыми элементами) ${{Q}_{{pv}}}$ и ${{W}_{{pv}}}$ на ${{\Upsilon }_{c}}$ могут быть выполнены с помощью транспонированной матрицы $R_{{pc}}^{T}$

Известно [17, 18], что спектр сеточной функции является периодическим и представляет собой бесконечный ряд сдвинутых копий спектра аппроксимируемой функции. Расстояние в частотной области между соседними копиями спектра равно частоте дискретизации ${{\omega }_{0}} = 2\pi {\text{/}}d$, где d – шаг сетки. Для исключения эффекта наложения частот и возможности восстановления периодической функции по ее сеточной аппроксимации будем полагать, что число гармоник в (7.5) и количество узлов сетки ${{\Upsilon }_{c}}$ удовлетворяют условиям

вытекающим из теоремы Котельникова–Шеннона [17, 18].

Введем матрицы Фурье ${{F}_{1}}$ и ${{F}_{2}}$ порядков ${{N}_{{c1}}}$ и ${{N}_{{c2}}}$ соответственно

Матрицы Фурье обратимы и при этом обратные матрицы имеют вид

где $F_{1}^{*},\;F_{2}^{*}$ – эрмитово сопряженные матрицы.

Используя матрицы F₁ и F₂, вычислим образы Фурье ${{\tilde {Q}}_{c}}$ и ${{\tilde {W}}_{c}}$ сеточных функций нормальных напряжений и перемещений

Для соответствующих матрицам ${{Q}_{c}}$ и ${{W}_{c}}$ векторов формулы (7.14) примут вид

где символом $ \otimes $ обозначено кронекерово произведение матриц ${{F}_{2}}$ и ${{F}_{1}}$ .

Сеточные функции ${{Q}_{c}}$ и ${{W}_{c}}$ являются вещественными, поэтому из свойств двумерного ДПФ [17, 18] следует, что элементы матриц ${{\tilde {Q}}_{c}}$ и ${{\tilde {W}}_{c}}$ удовлетворяют для ${{n}_{1}} = \overline {2,{{N}_{{c1}}}} $ и ${{n}_{2}} = \overline {2,{{N}_{{c2}}}} $ условиям

Для сеточных функций соотношение (7.3) примет вид

где символом $ \circ $ обозначено адамарово произведение матриц ${{\tilde {K}}_{s}}$ и ${{\tilde {Q}}_{c}}$.

Элементы матрицы ${{\tilde {K}}_{s}}$ с учетом (7.2) и (7.16) для индексов n₁ и n₂ из диапазонов $1 \leqslant {{n}_{1}} \leqslant {{N}_{{c1}}}{\text{/}}2 + 1$ и $1 \leqslant {{n}_{2}} \leqslant {{N}_{{c2}}}{\text{/}}2 + 1$ вычисляются по формуле

Остальные элементы матрицы ${{\tilde {K}}_{s}}$ определяются с помощью равенств

Несложно показать, что элемент

представляет собой коэффициент податливости рассматриваемого упругого слоя для случая приложения равномерной нормальной нагрузки по всей границе ${{\Gamma }_{1}}$.

Из (7.7)–(7.9), (7.14)–(7.15) и (7.17) следует, что сеточная аппроксимация ${{G}_{{s{v}}}}\,:\,{{Q}_{{p{v}}}}\, \mapsto \,{{W}_{{p{v}}}}$ оператора Пуанкаре–Стеклова ${{G}_{s}}$ имеет вид

где ${{\tilde {K}}_{{sv}}} = vec({{\tilde {K}}_{s}})$.

Учитывая, что все элементы ${{\tilde {K}}_{{s{v}}}}$ положительны, можно показать, что ${{G}_{{s{v}}}}$ является симметричной положительно-определенной квадратной матрицей порядка N. При численной реализации формировать матрицу ${{G}_{{s{v}}}}$ в явном виде не требуется, достаточно лишь программно реализовать вычисление произведений каждой из матриц в правой части (7.21) на векторы. Учитывая, что матрицы ${{R}_{{pc}}}$ и $R_{{pc}}^{T}$ содержат только N ненулевых элементов, причем равных 1, вычисление произведений этих матриц на векторы сводится к выполнению N операций присваивания. Умножение на вектор ${{\tilde {K}}_{{s{v}}}}$ требует выполнения N_c операций умножения. Для вычисления произведений матриц ${{F}_{2}} \otimes {{F}_{1}}$ и $({{F}_{2}} \otimes {{F}_{1}})_{{}}^{*}$ на векторы целесообразно использовать алгоритмы БПФ. Поэтому далее будем полагать ${{N}_{{c1}}} = {{2}^{{{{n}_{1}}}}},{{N}_{{c2}}} = {{2}^{{{{n}_{2}}}}},{{n}_{1}},{{n}_{2}} \in \mathbb{N}$. В этом случае общее число операций для выполнения одного преобразования имеет порядок $O({{N}_{c}}lo{{g}_{2}}{{N}_{c}})$ [17, 18].

Далее для аппроксимации задачи (6.5) применим гранично-элементный подход [35–37]. Для вычисления билинейной $g( \cdot , \cdot )$ и линейной $\phi ( \cdot )$ форм, определенных соответственно формулами (5.2) и (6.4), используем кубатурную формулу – прямое произведение квадратурных формул прямоугольников, узлы которой совпадают с узлами сетки ${{\Upsilon }_{p}}$. Тем самым определяются аппроксимирующие билинейная ${{g}_{h}}( \cdot , \cdot )$ и линейная ${{\phi }_{h}}( \cdot )$ формы в пространстве ${{\mathbb{R}}^{N}}$

где ${{\Phi }_{{v}}} = vec(\Phi ) \in {{\mathbb{R}}^{N}}$, $\Phi $ – сеточная аппроксимация на ${{\Upsilon }_{p}}$ функции $\Phi {\kern 1pt} ({{x}_{1}},{{x}_{1}})$, описывающей форму основания штампа.

При аппроксимации множества статически допустимых нормальных напряжений $\Sigma $, определенного формулой (6.1), применяется комбинированный подход. Для аппроксимации ограничений в виде неравенств используется коллокационный метод, а при аппроксимации ограничений в виде равенств – указанная выше кубатурная формула. В результате получается замкнутое выпуклое множество статически допустимых узловых нормальных напряжений

Коэффициенты ${{\lambda }_{i}}$, ${{\vartheta }_{i}}$ и ${{\theta }_{i}}$ вычисляются по формулам

где $x_{1}^{i},x_{2}^{i}$ – координаты i-го узла сетки ${{\Upsilon }_{p}}$.

В результате для задачи минимизации (6.5) получим сеточную аппроксимацию – задачу квадратичного программирования: найти сеточную функцию нормальных напряжений ${{Q}_{{p{v}}}} \in {{\mathbb{R}}^{N}}$ такую, что

Отметим, что размерность задачи (7.25) равна N – количеству узлов сетки ${{\Upsilon }_{p}}$ на Γ_p и не зависит от размера вычислительной области Γ_c, т.е. от выбора коэффициентов расширения вычислительной области ${{\chi }_{1}}$ и ${{\chi }_{2}}$.

8. Решение задачи квадратичного программирования. Одна из сложностей, возникающих при численном решении задачи (7.25), состоит в наличии ограничений в виде равенств, содержащих все переменные задачи. Аналогично [25] для упрощения вида ограничений введем линейное преобразование переменных $Z = CQ$ такое, что

Выбор индексов k, m и l в общем случае произволен при условии, что соответствующие компоненты вектора решения ${{Q}_{{p{v}}}}$ отличны от нуля, т.е. выполняются неравенства ${{q}_{k}} < 0,$ ${{q}_{m}} < 0$ и ${{q}_{l}} < 0$. При решении конкретной задачи такие индексы (номера узлов) могут быть выбраны из априорных соображений. Корректность сделанного выбора проверяется путем апостериорного анализа полученного решения.

Нетрудно видеть, что матрица C является невырожденной. Выражения для элементов матриц ${{C}^{{ - 1}}}$ и ${{C}^{{ - T}}}$ несложно получить в явном виде, однако формировать матрицы в явном виде не требуется, достаточно лишь программно реализовать вычисление произведений этих матриц на векторы.

В результате преобразования переменных (8.1)–(8.2) получим эквивалентную (7.25) задачу квадратичного программирования: найти элемент ${{Z}_{p}} \in {{V}_{z}}$ такой, что

где

Задача (8.3) является задачей квадратичного программирования с простыми ограничениями. Учитывая, что решение задачи (7.25) не зависит от значения множителя ${{d}_{1}}{{d}_{2}}$ – площади граничного элемента, в (8.3) этот множитель опущен.

Из (8.4) следует, что матрицы ${{G}_{{s{v}}}}$ и A_z подобны, поэтому матрица Гессе A_z минимизируемой функции ${{J}_{z}}(Z)$ также является симметричной положительно-определенной матрицей.

Для численного решения задачи квадратичного программирования (8.3) в настоящей работе используется вариант МСГ, подробно рассмотренный в [25] при решении задач дискретного контакта для упругой полуплоскости.

9. Численные результаты. Разработанный вычислительный алгоритм реализован на языке FORTRAN с применением программного пакета для разработчиков NVIDIA HPC SDK. Для выполнения двумерных ДПФ использовалась библиотека cuFFT, позволяющая с помощью технологии CUDA производить вычисления на графических процессорах.

Для верификации алгоритма и разработанного программного обеспечения проведено сравнение численных решений ряда задач дискретного контакта, в частности, рассмотренных ниже задач для штампов с регулярным поверхностным рельефом с решениями, полученными с помощью известного алгоритма [24]. Учитывая возможности последнего, при постановке контактных задач использовалось лишь первое из уравнений равновесия штампа (2.5) и полагалось ${{\varphi }_{1}} = {{\varphi }_{1}} = 0$. Проведенные расчеты показали, что среднеквадратичные относительные расхождения решений (сеточных функций нормальных напряжений) не превышают 10^–5.

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

где ${{b}_{1}},{{b}_{2}}$ – размеры области контакта ${{\Gamma }_{p}}$, ${{\alpha }_{1}},{{\alpha }_{2}} \geqslant 0$ – безразмерные параметры.

Несложно показать аналитически, что в случае нагружения такого штампа с эксцентриситетом $({{e}_{1}}{{b}_{1}},{{e}_{2}}{{b}_{2}})$ относительно точки начального контакта ${{x}_{1}} = {{x}_{2}} = 0$ распределение нормальных напряжений будет соответствовать распределению при отсутствии эксцентриситета, смещенному вдоль осей координат $O{{x}_{1}}$ и $O{{x}_{2}}$ соответственно на ${{e}_{1}}{{b}_{1}}$ и ${{e}_{2}}{{b}_{2}}$, осадка штампа ${{\delta }_{3}}$ уменьшится на ${{\alpha }_{1}}{{b}_{1}}e_{1}^{2} + {{\alpha }_{2}}{{b}_{2}}e_{2}^{2}$, а сам штамп повернется на углы ${{\varphi }_{1}} = - 2{{\alpha }_{2}}{{e}_{2}}$ и ${{\varphi }_{2}} = 2{{\alpha }_{1}}{{e}_{1}}$. Проведенные расчеты подтвердили эти выводы.

Далее приведен анализ результатов решения задач дискретного контакта, в которых область фактического контакта состояла из совокупности отдельных пятен контакта. Рассматривались задачи для гладких прямоугольных в плане штампов с регулярным поверхностным рельефом, состоящим из $K = {{K}_{1}} \times {{K}_{2}}$ микровыступов. Номинальная область контакта полагалась равной ${{\Gamma }_{p}} = \{ 0 \leqslant {{x}_{1}} \leqslant {{b}_{1}},0 \leqslant {{x}_{2}} \leqslant {{b}_{2}},{{x}_{3}} = h\} $. Форма основания штампов задавалась функцией

где $\Phi {\kern 1pt} {\text{*}}({{x}_{1}},{{x}_{2}})$ – выпуклая функция, определяющая макроформу штампа; ${{\Psi }_{1}}({{\xi }_{1}})$, ${{\Psi }_{2}}({{\xi }_{2}})$ – строго выпуклые функции, характеризующие форму микровыступов; ${{\xi }_{1}} = \left\{ {{{K}_{1}}{{x}_{1}}{\text{/}}{{b}_{1}}} \right\}$, ${{\xi }_{2}} = \left\{ {{{K}_{2}}{{x}_{2}}{\text{/}}{{b}_{2}}} \right\}$ – “быстрые” координаты, $\left\{ \cdot \right\}$ – дробная часть числа.

Для заданной тройки функций $\Phi {\text{*}}({{x}_{1}}, {{x}_{2}})$, ${{\Psi }_{1}}({{\xi }_{1}})$ и ${{\Psi }_{2}}({{\xi }_{2}})$ формула (9.2) определяет двухпараметрическое семейство штампов $\Xi ({{{\rm K}}_{1}},{{{\rm K}}_{2}})$, в качестве параметров которого выступают числа микровыступов ${{K}_{1}}$ и ${{K}_{2}}$ вдоль сторон штампа. Штампы, принадлежащие к одному семейству, имеют одинаковую макроформу $\Phi {\kern 1pt} {\text{*}}({{x}_{1}},{{x}_{2}})$, а поперечные сечения их микровыступов, параллельные соответствующим координатным плоскостям, являются подобными.

Расчеты выполнялись для следующего семейства жестких штампов:

где ${{\alpha }_{1}},{\kern 1pt} \,{{\alpha }_{2}} \geqslant 0,$ ${{\beta }_{1}},{\kern 1pt} \,{{\beta }_{2}} > 0,$ ${{m}_{1}},\,{{m}_{2}},{{n}_{1}},\,{{n}_{2}} > 1$ – безразмерные параметры.

Компоненты главного вектора и главного момента внешних сил, приложенных к штампу, задавались в виде:

где $E{\kern 1pt} {\text{*}} = E{\text{/(}}1 - {{\nu }^{2}}{\text{)}}$ – приведенный модуль упругости; f > 0 – безразмерный параметр; ${{e}_{1}},{{e}_{2}}$ – безразмерные параметры, характеризующие эксцентриситет равнодействующей внешней нагрузки относительно центра приведения ${{x}^{c}} = (0.5\,{{b}_{1}},0.5\,{{b}_{2}},h)$. При проведении расчетов значения параметров внешней нагрузки  f, ${{e}_{1}}$ и ${{e}_{2}}$ выбирались таким образом, чтобы решение контактной задачи существовало и пятна контакта отдельных микровыступов оставались изолированными друг от друга.

Введем на множестве микровыступов сеточные координаты $(i,j),$ $i = \overline {1,{{K}_{1}}} ,$ $j = \overline {1,{{K}_{2}}} ,$ и обозначим через ${{\Gamma }_{{ij}}}$ часть ${{\Gamma }_{p}}$, соответствующую микровыступу с координатами $(i,j)$. Отметим, что каждый штамп рассматриваемого семейства $\Xi ({{{\rm K}}_{1}},{{{\rm K}}_{2}})$ имеет свою собственную систему сеточных координат, определяемую параметрами ${{K}_{1}}$ и ${{K}_{2}}$. Далее для каждого штампа введем следующие пары сеточные функций:

– значений $P = [{{p}_{{ij}}}]$ и координат ${{X}_{p}} = [{{({{x}_{1}},{{x}_{2}})}_{{ij}}}]$ максимумов контактного давления на микровыступах

– нормализованных контактных усилий $R = [{{r}_{{ij}}}]$ и координат ${{X}_{r}} = [{{({{x}_{1}},{{x}_{2}})}_{{ij}}}]$ точек их приложения на микровыступах – относительных величин площадей $S = [{{s}_{{ij}}}]$ и координат ${{X}_{s}} = [{{({{x}_{1}},{{x}_{2}})}_{{ij}}}]$ центров тяжести пятен фактического контакта на микровыступах где $[ \cdot ]$ – скобка Айверсона (функция равная 1 для истинного аргумента и равная 0 в противном случае).

Для каждой пары функций $\{ P,{{X}_{p}}\} $, $\{ R,{{X}_{r}}\} $ и $\{ S,{{X}_{s}}\} $ с помощью некоторого оператора интерполяции построим на Γ_p непрерывные функции $p({{x}_{1}},{{x}_{2}})$, $r({{x}_{1}},{{x}_{2}})$ и $s({{x}_{1}},{{x}_{2}})$, называемые далее соответственно огибающей контактного давления, огибающей нормализованных контактных усилий и огибающей относительных величин фактических площадей контакта микровыступов.

В результате обработки результатов многочисленных вычислительных экспериментов установлена следующая закономерность: если условия нагружения штампов таковы, что пятна контакта отдельных микровыступов остаются изолированными друг от друга, то для двухпараметрического семейства штампов $\Xi ({{{\rm K}}_{1}},{{{\rm K}}_{2}})$ существуют единая огибающая контактного давления $\tilde {p}({{x}_{1}},{{x}_{2}})$, единая огибающая нормализованных контактных усилий $\tilde {r}({{x}_{1}},{{x}_{2}})$ и единая огибающая относительных величин фактических площадей контакта $\tilde {s}({{x}_{1}},{{x}_{2}})$ микровыступов, форма и положение которых зависят от параметров f, ${{e}_{1}}$ и ${{e}_{2}}$ внешней нагрузки, а также отношения размеров области контакта Γ_p к толщине упругого слоя h. В качестве $\tilde {p}({{x}_{1}},{{x}_{2}})$, $\tilde {r}({{x}_{1}},{{x}_{2}})$ и $\tilde {s}({{x}_{1}},{{x}_{2}})$ можно выбирать соответственно огибающие $p({{x}_{1}},{{x}_{2}})$, $r({{x}_{1}},{{x}_{2}})$ и $s({{x}_{1}},{{x}_{2}})$ для некоторого “базового” штампа рассматриваемого семейства $\Xi ({{{\rm K}}_{1}},{{{\rm K}}_{2}})$, например, штампа с наибольшим числом микровыступов $K = {{K}_{1}} \cdot {{K}_{2}}$ в проведенной серии расчетов. Альтернативный подход к построению единых огибающих состоит в применении некоторой процедуры осреднения данных по всем штампам рассматриваемого семейства $\Xi ({{{\rm K}}_{1}},{{{\rm K}}_{2}})$.

В качестве примера в табл. 1 приведены среднеквадратичные относительные отклонения

для семейства штампов $\Xi ({{{\rm K}}_{1}},{{{\rm K}}_{2}})$ при следующих значениях параметров: ${{\alpha }_{1}}\, = \,{{\alpha }_{2}}$ = 3 × 10^–5, ${{\beta }_{1}} = {{\beta }_{2}} = 1{{0}^{{ - 4}}}$, ${{m}_{1}} = {{m}_{2}} = 2$, ${{n}_{1}} = {{n}_{2}} = 4$, ${{b}_{1}} = {{b}_{2}} = b$. В качестве базового выбирался штамп со 128 × 128 микровыступами. Параметры внешней нагрузки в (9.5) полагались следующими:  f = 2 × 10^–5, ${{e}_{1}} = 0.05$ и ${{e}_{2}} = 0$. Толщина упругого слоя принималась равной $h = b{\text{/}}5$.

Сеточные нормы ${{\left\| \cdot \right\|}_{2}}$ в (9.9) вычислялись соответственно в узлах сеток ${{X}_{p}}$, ${{X}_{r}}$ и ${{X}_{s}}$ для каждого штампа рассматриваемого семейства. Интерполяция функций $\tilde {p}({{x}_{1}},{{x}_{2}})$, $\tilde {r}({{x}_{1}},{{x}_{2}})$ и $\tilde {s}({{x}_{1}},{{x}_{2}})$ производилась с помощью метода обратно взвешенных расстояний.

В табл. 2 приведены аналогичные данные для семейства штампов $\Xi ({{{\rm K}}_{1}},{{{\rm K}}_{2}})$ при следующих значениях параметров: ${{m}_{1}} = 2$, ${{m}_{2}} = 4$, ${{n}_{1}} = 4$, ${{n}_{2}} = 2$, ${{e}_{1}} = {{e}_{2}} = 0.05$, $h = b{\text{/}}2$. Значения остальных параметров принимались такими же, как и для табл. 1.

Отметим, что при проведении расчетов необходимое количество $N = {{N}_{{p1}}} \cdot {{N}_{{p2}}}$ граничных элементов определялось путем сравнения решений, полученных на вложенных сетках при их двукратном последовательном измельчении. При решении рассмотренных выше задач для штампов с регулярным поверхностным рельефом каждый микровыступ ${{\Gamma }_{{ij}}}$ аппроксимировался сеткой из 64 × 64 граничных элементов. Наибольшее количество элементов сетки для базового штампа со 128 × 128 микровыступами составляло 2²⁶ = 67 108 864 граничных элемента.

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

Оператор Стеклова–Пуанкаре для упругого слоя строится с помощью двумерного интегрального преобразования Фурье, поэтому при аппроксимации этого оператора используется двумерное ДПФ. Для численной реализации прямого и обратного ДПФ применяются алгоритмы БПФ.

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

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