A unilateral discrete contact problem for a stratified elastic strip

Full Text

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

При исследовании локального контактного взаимодействия тел с покрытиями в качестве расчетной схемы, как правило, выбирается упругая полоса или слой, сцепленные с основанием. Учет поверхностного микрорельефа контактирующих тел в виде волнистости или шероховатости приводит к постановке задач дискретного (множественного) контакта [2$–$6]. Подробный обзор современного состояния исследований в области механики дискретного контакта, включая основные подходы к постановке задач, методы аналитического и численного решения, конкретные результаты и области их практического использования, приведен в статье [7].

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

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

2. Постановка задачи. Пусть невесомая стратифицированная упругая полоса в прямоугольной системе координат $O{x}_1{x}_2$ занимает область $\Omega =\{x=(x_1,x_2)\in {ℝ}^2:\left|x_1\right|\le \infty ,0\le x_2\le h\}$ и состоит из произвольного числа $N$ изотропных упругих слоев ${\Omega }_n$, границы раздела которых параллельны оси $O{x}_1$. Слои пронумерованы как $\mathrm{1,}\dots ,N$ в порядке возрастания координаты $x_2$. Границу полосы $x_2\equiv 0$ обозначим ${\Gamma }_0$, границу раздела слоев ${\Omega }_n$ и ${\Omega }_{n+1}$ $–$ через ${\Gamma }_n$, а границу полосы $x_2\equiv h$ $–$ через ${\Gamma }_N$. Параметры Ламе материала полосы являются произвольными ограниченными функциями координаты $x_2$: $\lambda =\lambda (x_2)$ и $\mu =\mu (x_2)$, имеющими разрывы первого рода на границах ${\Gamma }_n$ раздела слоев. Из физических соображений следует, что существуют постоянные ${\lambda }_0>0$ и ${\mu }_0>0$, такие, что

$\lambda (x_2)\ge {\lambda }_0$, $\mu (x_2)\ge {\mu }_0$, $0\le x_2\le h$ (2.1)

Далее под $u^{(n)}(x),{\epsilon }^{(n)}(x),{\sigma }^{(n)}(x)$ будем понимать соответственно вектор перемещений и тензоры деформаций и напряжений в точке $x\in {\Omega }_n$.

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

Предполагается, что упругая полоса находится в условиях плоской деформации, деформации малы, а массовые силы и напряжения в недеформированном состоянии отсутствуют. Напряженно-деформированное состояние слоев ${\Omega }_n,$ $n=\overline{\mathrm{1,}N},$ описывается системой уравнений:

$\epsilon =defu$, $\sigma =S:\epsilon$, $div\sigma =0$ в ${\Omega }_n$, (2.2)

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

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

$u_1=u_2=0$ на ${\Gamma }_0$ (2.3)

На границах раздела слоев также задаются условия полного сцепления:

$u_1^{(n)}=u_1^{(n+1)},$ $u_2^{(n)}=u_2^{(n+1)},$ ${\sigma }_{12}^{(n)}={\sigma }_{12}^{(n+1)},$ ${\sigma }_{22}^{(n)}={\sigma }_{22}^{(n+1)}$ на ${\Gamma }_n$ (2.4)

В упругую полосу вдавливается гладкий жесткий штамп, основание которого имеет поверхностный микрорельеф. Часть границы ${\Gamma }_N$, по которой возможен контакт полосы со штампом, обозначается ${\Gamma }_p$. Положение и предельные размеры ${\Gamma }_p$, т.е. номинальная область контакта, задаются априори, исходя из геометрических соображений. Предполагается, что часть границы ${\Gamma }_p$ является односвязной и конечной. При вдавливании штампа с поверхностным микрорельефом номинальная область контакта ${\Gamma }_p$ включает в себя множество отдельных пятен фактического контакта, положение и размеры которых заранее неизвестны.

Форма основания штампа и его поверхностный микрорельеф описываются функцией $\Phi \text{}(x_1)$, значение которой в точке $x\in {\Gamma }_p$ равно расстоянию от этой точки до поверхности штампа, измеренному вдоль направления внешней нормали к границе ${\Gamma }_p$. Расстояние $\Phi (x_1)$ отсчитывается по отношению к недеформированному состоянию полосы. Для определенности будем полагать $\underset{x\in {\Gamma }_p}{min}\Phi (x_1)=0$. Для штампа с поверхностным микрорельефом функция $\Phi (x_1)$ является мультимодальной (многоэкстремальной). Положение штампа определяется вектором перемещений $\delta =({\delta }_1,{\delta }_2)$ и углом поворота ${\phi }_3$ штампа как жесткого целого. Перемещения и углы поворота штампа предполагаются малыми. Главный вектор $F=(F_1,F_2)$ и главный момент $M_3$ внешних сил, приложенных к штампу, считаются заданными. В качестве центра приведения выбирается точка $x^c=(x_1^c,x_2^c)$. Далее рассматривается задача нормального контакта упругой полосы со штампом, поэтому будем полагать

${\delta }_1=0,$ $F_1=0,$ $\infty <F_2<0,$ $\left|M_3\right|<\infty$

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

$u_2\le \Phi +{\delta }_2+(x_1-x_1^c){\phi }_3$, ${\sigma }_{22}\le 0$, ${\sigma }_{12}=0$

${\sigma }_{22}[u_2-\Phi -{\delta }_2-(x_1-x_1^c){\phi }_3]=0$ на ${\Gamma }_p$ (2.5)

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

${\sigma }_{21}={\sigma }_{22}=0$ на ${\Gamma }_N\backslash {\Gamma }_p$ (2.6)

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

${\displaystyle \underset{{\Gamma }_p}{\int }{\sigma }_{22}}d{\Gamma }_p=F_2$, ${\displaystyle \underset{{\Gamma }_p}{\int }(x_1-x_1^c){\sigma }_{22}}d{\Gamma }_p=M_3$  (2.7)

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

Для существования решения рассматриваемой контактной задачи далее будем предполагать, что внешние силы и моменты, приложенные к жесткому штампу, согласованы между собой таким образом, что существует распределение нормальных напряжений ${\sigma }_{22}\le 0$ на ${\Gamma }_p$, удовлетворяющее уравнениям равновесия штампа (2.7).

Для выделения класса единственности решения в рассматриваемой контактной задаче используется условие конечности потенциальной энергии деформации упругой полосы [13]:

${\displaystyle \underset{\Omega }{\int }\sigma :\epsilon }d\Omega <\infty$ (2.8)

Задача (в дифференциальной постановке) состоит в определении полей перемещений $u^{(n)}(x)$, деформаций ${\epsilon }^{(n)}(x)$ и напряжений ${\sigma }^{(n)}(x)$, $n=\overline{\mathrm{1,}N},$ удовлетворяющих уравнениям (2.2), граничным условиям (2.3)$–$(2.6), уравнениям равновесия штампа (2.7) и условию (2.8). Также необходимо найти смещение ${\delta }_2$ и поворот ${\phi }_3$ штампа. Подчеркнем, что в рассматриваемой задаче одностороннего дискретного контакта априори задается лишь номинальная область контакта ${\Gamma }_p$, положение и размеры пятен фактического контакта заранее неизвестны и подлежат определению в процессе решения задачи.

3. Интегральное представление решения. Рассмотрим вспомогательную краевую задачу для стратифицированной упругой полосы $\Omega$: найти поля перемещений $u^{(n)}$, $n=\overline{\mathrm{1,}N},$ удовлетворяющие уравнениям (2.2), граничным условиям (2.3), (2.4) и (2.6), условию (2.8), а также граничным условиям

${\sigma }_{21}=0$, ${\sigma }_{22}=q$ на ${\Gamma }_p$

Из результатов [13$–$16] следует, что если выполняются условия (2.1), то для любого $q\in L_2({\Gamma }_p)$ данная задача имеет единственное решение, которое может быть представлено в виде

$u(x)={\displaystyle \underset{{\Gamma }_p}{\int }g(x,{\xi }_1)q({\xi }_1)d\text{}{\xi }_1}$; $x\in {\Omega }_n$, $n=\overline{\mathrm{1,}N}$, (3.1)

где $g(x,{\xi }_1)$ $–$ ядро интегрального представления, которое может быть интерпретировано как поле перемещений $–$ решение краевой задачи о действии нормальной сосредоточенной единичной силы, приложенной в точке $x=({\xi }_1,h)$ границы ${\Gamma }_p$ стратифицированной упругой полосы $\Omega$.

Решение (3.2) вспомогательной краевой задачи удовлетворяет всем условиям сформулированной в разд. 2 задачи одностороннего дискретного контакта, кроме условий одностороннего гладкого контакта (2.5) и условий равновесия штампа (2.7). Поэтому использование интегрального представления (3.1) позволяет свести решение рассматриваемой контактной задачи к нахождению на ${\Gamma }_p$ нормальных напряжений $q={\sigma }_{22}\in L_2({\Gamma }_p)$, удовлетворяющих следующей системе уравнений и неравенств:

$q\le 0$ на ${\Gamma }_p$ (3.2)

$u_2(q)\le \Phi +{\delta }_2+(x_1-x_1^c){\phi }_3$ на ${\Gamma }_p$ (3.3)

$q[u_2(q)-\Phi -{\delta }_2-(x_1-x_1^c){\phi }_3]=0$ на ${\Gamma }_p$ (3.4)

${\displaystyle \underset{{\Gamma }_p}{\int }q}d{\Gamma }_p=F_2$, ${\displaystyle \underset{{\Gamma }_p}{\int }(x_1-x_1^c)q}d{\Gamma }_p=M_3$, (3.5)

где $u_2(q)$ $–$ нормальные перемещения на ${\Gamma }_p$, соответствующие нормальным напряжениям $q$.

4. Оператор Пуанкаре$–$Стеклова. Для решения системы (3.2)$–$(3.5) построим ОПС $S_{ps}:q\mapsto w$, отображающий посредством решения (3.1) нормальные напряжения $q(x_1)\equiv {\sigma }_{22}(x_1,h)$ на части ${\Gamma }_p$ границы упругой полосы в нормальные перемещения $w(x_1)\equiv u_2(x_1,h)$ на ${\Gamma }_p$.

Соответствующее решению (3.1) представление ОПС $S_{ps}:q\mapsto w$ для $q\in L_2({\Gamma }_p)$ имеет вид:

$w(x_1)={\displaystyle \underset{{\Gamma }_p}{\int }{g}_{ps}(x_1-{\xi }_1)q({\xi }_1)d\text{}{\xi }_1}$ (4.1)

Выражение в явном виде ядра $g_{ps}(\cdot )$ интегрального представления (4.1) для произвольной стратифицированной упругой полосы неизвестно. Правая часть (4.1) является интегральным оператором типа свертки, поэтому с помощью интегрального преобразования Фурье по координате $x_1$ можно получить алгебраическое соотношение, связывающее трансформанты нормальных перемещений $\tilde{w}(\alpha )$ и нормальных напряжений $\tilde{q}(\alpha )$

$\tilde{w}(\alpha )={\tilde{g}}_{ps}(\alpha )\tilde{q}(\alpha )$, (4.2)

где $\alpha$ $–$ параметр преобразования Фурье, ${\tilde{g}}_{ps}(\alpha )$ $–$ трансформанта ядра интегрального представления (4.1). Таким образом действие ОПС $S_{ps}$ сводится к выполнению прямого и обратного преобразований Фурье и перемножению трансформант.

5. Передаточная функция оператора Пуанкаре$–$Стеклова. Функцию ${\tilde{g}}_{ps}(\alpha )$, следуя [17], будем называть передаточной функцией. В случае однородной полосы выражение для передаточной функции может быть получено аналитически [13]. Для функционально-градиентной полосы эту функцию удается построить с помощью численно-аналитической методики при специальной зависимости ее упругих свойств по толщине, в частности степенной или экспоненциальной зависимости [4]. В случае произвольного закона изменения упругих свойств по толщине используются приближенные подходы, основанные как на прямом численном интегрировании краевых задач для систем дифференциальных уравнений по поперечной координате [15, 17], так и на замене неоднородной полосы многослойной с кусочно-постоянной зависимостью упругих модулей от координаты и сведении задачи к решению систем функциональных уравнений [16].

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

В настоящей работе для построения передаточной функции применялся альтернативный подход [12], опирающийся на использование вариационной формулировки краевой задачи для трансформант перемещений. Аппроксимация вариационных уравнений производилась методом конечных элементов. Для численного решения задачи использован нестационарный итерационный метод, на каждом шаге которого методом прогонки решались две системы линейных алгебраических уравнений с трехдиагональными матрицами. Разработан алгоритм выбора последовательности параметров итерационного метода, обеспечивающей его сходимость для любых значений параметра преобразования Фурье $\alpha$.

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

При решении задач дискретного контакта требуется вычислять значения передаточной функции ОПС в достаточно широком диапазоне изменения параметра преобразования Фурье $\alpha$ [8]. С целью уменьшения вычислительных затрат в настоящей работе для больших значений параметра $\alpha$ применялся алгоритм, использующий полученное в [19] трехчленное асимптотическое разложение передаточной функции и построенные на его основе аппроксимации Паде.

6. Вариационная формулировка задачи. Для решения системы (3.2)$–$(3.5), содержащей уравнения и неравенства, используется вариационный подход [9]. Введем необходимые для построения вариационных формулировок задачи функциональные пространства $–$ гильбертов триплет $H^{1/2}({\Gamma }_p)\subset L_2({\Gamma }_p)\subset {\tilde{H}}^{-1/2}({\Gamma }_p)$ [8]. Учитывая плотность вложения $L_2({\Gamma }_p)\subset {\tilde{H}}^{-1/2}({\Gamma }_p)$, ОПС $S_{ps}$, определенный формулой (4.1), можно продолжить по непрерывности на ${\tilde{H}}^{-1/2}({\Gamma }_p)$ и аналогично [8] рассматривать как оператор, действующий из ${\tilde{H}}^{-1/2}({\Gamma }_p)$ в $H^{1/2}({\Gamma }_p)$.

Обозначим через $⟨\cdot ,\cdot ⟩$ отношение двойственности ${}_{{\tilde{H}}^{-1/2}({\Gamma }_p)}⟨\cdot ,\cdot ⟩{}_{H^{1/2}({\Gamma }_p)}$, порожденное продолжением скалярного произведения в $L_2({\Gamma }_p)$. Введем на ${\tilde{H}}^{-1/2}({\Gamma }_p)$ непрерывные билинейную и линейную формы

$a(p,q)=⟨p,G_{ps}q⟩$, $b(p)=⟨p,\Phi ⟩$,

предполагая, что $\Phi \in H^{1/2}({\Gamma }_p)$. Аналогично [8] можно показать, что если выполняются условия (2.1), то билинейная форма $a(\cdot ,\cdot )$ является симметричной и коэрцитивной на ${\tilde{H}}^{-1/2}({\Gamma }_p)$.

Образуем множество статически допустимых нормальных напряжений на ${\Gamma }_p$

$V=\left\{p\in {\tilde{H}}^{-1/2}({\Gamma }_p):p\le 0;⟨p,1⟩=F_2;⟨p,x_1-x_1^c⟩=M_3\right\}$ (6.1)

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

Аналогично [8] можно показать, что решение $q\in L_2({\Gamma }_p)$ системы (3.2)$–$(3.5) удовлетворяет граничному вариационному неравенству: найти $q\in V$ такой, что

$a(p-q,q)-b(p-q)\ge 0\forall p\in V$ (6.2)

Отметим, что неравенство (6.2) не содержит неизвестных смещения ${\delta }_2$ и угла поворота ${\phi }_3$ жесткого штампа. Как показано в [8], это является следствием того, что элементы множества статически допустимых нормальных напряжений $V$ удовлетворяют уравнениям равновесия штампа (3.5).

Используя известные приемы [9, 20, 21], можно доказать обратное утверждение: решение $q$ вариационного неравенства (6.2) удовлетворяет системе уравнений и неравенств (3.2)$–$(3.5) по крайней мере в обобщенном смысле, т.е. условия (3.2) и (3.4) выполняются в смысле пространства ${\tilde{H}}^{-1/2}({\Gamma }_p)$, условие (3.3) выполняется в смысле пространства $H^{1/2}({\Gamma }_p)$, а условия (3.5) $–$ после продолжения по непрерывности на ${\tilde{H}}^{-1/2}({\Gamma }_p)$:

$⟨q,1⟩=F_2$, $⟨q,x_1-x_1^c⟩=M_3$

Учитывая, что билинейная форма $a(\cdot ,\cdot )$ является положительно определенной, а множество $V$ $–$ непустым замкнутым выпуклым множеством в ${\tilde{H}}^{-1/2}({\Gamma }_p)$, вариационное неравенство (6.2) эквивалентно задаче минимизации граничного функционала: найти $q\in V$ такой, что

$J(q)=\underset{p\in V}{inf}\left\{J(p)=a(p,p)/2-b(p)\right\}$ (6.3)

Кроме того, решение вариационного неравенства (6.2) и задачи минимизации (6.3) существует и единственно [9, 21].

7. Конечномерная аппроксимация задачи. Выбор метода аппроксимации задачи (6.3) определяется, в первую очередь, возможностью построения эффективной с вычислительной точки зрения аппроксимации ОПС (4.1). Из (4.2) следует, что действие ОПС сводится к выполнению пары (прямого и обратного) преобразований Фурье и перемножению трансформант. Использование преобразования Фурье наиболее эффективно в случае периодических функций благодаря дискретности их спектра и, как следствие, переходе от непрерывного преобразования к дискретному. Учитывая, что для стратифицированной полосы передаточная функция ОПС строится численно, основная идея используемого подхода состоит в аппроксимации искомых нормальных напряжений периодическими сеточными функциями и применении алгоритмов быстрого преобразования Фурье (БПФ). Для уменьшения возникающей при этом ошибки периодичности [22, 23] вводится расширенная вычислительная область ${\Gamma }_c\subset {\Gamma }_N$ такая, что

${\Gamma }_p⋐{\Gamma }_c,D=diam{\Gamma }_c=k_{c}diam{\Gamma }_p;k_c\in \mathbb{N},k_c>1$, (7.1)

где $k_c$ $–$ коэффициент расширения вычислительной области.

Построим на ${\Gamma }_c$ регулярную (равномерную) сетку $T_c$, состоящую из $M$ одноузловых граничных элементов нулевого порядка размером ${\delta }_e=D/M$. Узлы $T_c$ обозначим $t_m,m=\overline{1,M}$. Часть $T_c$, расположенную на ${\Gamma }_p$ и содержащую $L$ элементов, обозначим $T_p$. Из (7.1) следует, что $L=M/k_c$.

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

Введем на $T_c$ сеточные функции нормальных напряжений и перемещений

$Q_c=(q(t_1),q(t_2),\dots ,q(t_M))$, $W_c=(w(t_1),w(t_2),\dots ,w(t_M))$,

а также редукции этих функций на $T_p$

$Q_p=R{Q}_c$, $W_p=R{W}_c$,

где $R$ $–$ прямоугольная матрица размеров $L\times M$, в каждой строке которой имеется равно один ненулевой элемент, равный единице.

Обратные операции тривиального продолжения (дополнения нулевыми элементами) $Q_p$ и $W_p$ на $T_c$ могут быть выполнены с помощью транспонированной матрицы $R^T$

$Q_c=R^T{Q}_p$, $W_c=R^T{W}_p$

Введем матрицу Фурье порядка $M$

$F=\left[f_{mk}\right],f_{mk}={\omega }^{(m-1)(k-1)},\omega =exp\left(-2\pi i/M\right)$

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

$F^{-1}=F^{*}/M$,

где $F^{*}$ $–$ эрмитово сопряженная матрица.

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

${\tilde{Q}}_c=({\tilde{q}}_0,{\tilde{q}}_1,\dots ,{\tilde{q}}_{M-1})=F{Q}_c=F{R}^T{Q}_p$ (7.2)

${\tilde{W}}_c=({\tilde{w}}_0,{\tilde{w}}_1,\dots ,{\tilde{w}}_{M-1})=F{W}_c=F{R}^T{W}_p$ (7.3)

Сеточные функции $Q_c$ и $W_c$ являются вещественными, поэтому из свойств ДПФ следует, что компоненты ${\tilde{Q}}_c$ и ${\tilde{W}}_c$ удовлетворяют условиям

${\tilde{q}}_{M-m}={\tilde{q}}_m^{*},{\tilde{w}}_{M-m}={\tilde{w}}_m^{*};m=1,2,\dots ,M-1$ (7.4)

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

${\tilde{W}}_c=\tilde{G}{\tilde{Q}}_c;\tilde{G}=diag({\tilde{g}}_0,{\tilde{g}}_1,\dots ,{\tilde{g}}_{M-1})$,  (7.5)

где с учетом (7.4) для четных $M$

${\tilde{g}}_0={\tilde{g}}_{ps}(0)={\displaystyle \underset{0}{\overset{h}{\int }}\frac{d{x}_2}{\lambda (x_2)+2\mu (x_2)}}$, ${\tilde{g}}_m={\tilde{g}}_{M-m}={\tilde{g}}_{ps}\left(2\pi m/D\right);m=1,2,\dots ,M/2$

Несложно показать, что ${\tilde{g}}_0$ представляет собой коэффициент податливости рассматриваемой стратифицированной упругой полосы для случая приложения равномерной нормальной нагрузки по всей границе ${\Gamma }_N$.

Из (7.2), (7.3) и (7.5) следует, что сеточная аппроксимация интегрального представления (4.1), определяющего ОПС $S_{ps}$, имеет вид

$W_p=S{Q}_p$, (7.6)

где

$S=R{F}^{*}\tilde{G}F{R}^T/M$ (7.7)

Матрица $S$ является квадратной матрицей порядка $L$. Для вычисления произведения этой матрицы на вектор формировать ее в явном виде не требуется, достаточно программно реализовать вычисление произведений каждой из матриц в правой части (7.7) на векторы. Учитывая, что матрицы $R$ и $R^T$ содержат только $L$ ненулевых элементов равных $1$, вычисление произведений этих матриц на векторы сводится к выполнению $L$ операций присваивания. Умножение диагональной матрицы $\tilde{G}$ на вектор требует выполнения $M$ операций умножения. Для вычисления произведений матриц $F$ и $F^{*}$ на векторы целесообразно использовать алгоритмы БПФ. Поэтому далее будем полагать $M=2^m,m\in \mathbb{N}$. В этом случае число операций для выполнения одного преобразования имеет порядок $O(Mm)=O(Mlo{g}_{2}M)$.

Аппроксимация задачи минимизации (6.3) производится с помощью гранично-элементного подхода. Для вычисления билинейной $a(\cdot ,\cdot )$ и линейной $b(\cdot )$ форм используется квадратурная формула прямоугольников, узлы которой совпадают с узлами сетки $T_p$. Тем самым определяются аппроксимирующие билинейная $a_1(\cdot ,\cdot )$ и линейная $b_1(\cdot )$ формы в пространстве ${ℝ}^L$

$a_1(P_p,Q_p)=P_p^{T}S{Q}_p{\delta }_e,b_1(P_p)={\Phi }^T{P}_p{\delta }_e$,

где $\Phi ={({\varphi }_1,{\varphi }_2,\dots ,{\varphi }_L)}^T$ $–$ сеточная аппроксимация на $T_p$ функции $\Phi \text{}(x_1)$, описывающей форму основания штампа.

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

$V_1=\left\{Q\in {ℝ}^L:q_i\le 0,i=\overline{1,L};{\displaystyle \sum _{i=1}^L{\theta }_i{q}_i}=F_2;{\displaystyle \sum _{i=1}^L{\vartheta }_i{q}_i}=M_3\right\}$

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

${\theta }_i={\delta }_e,{\vartheta }_i={\delta }_e(x_1^i-x_1^c);i=\overline{1,L}$,

где $x_1^i$ $–$ координата узла $t_i$ сетки $T_p$.

В результате для задачи минимизации (6.3) получим сеточную аппроксимацию $–$ задачу квадратичного программирования: найти сеточную функцию нормальных напряжений $Q_p\in {ℝ}^L$ такую, что

$J_1(Q_p)=\underset{Q\in V}{min}\left\{J_1(Q)=Q^{T}SQ{\delta }_e/2-{\Phi }^{T}Q{\delta }_e\right\}$ (7.8)

Отметим, что размерность задачи (7.8) равна $L$ $–$ количеству узлов сетки на ${\Gamma }_p$ и не зависит от размера вычислительной области ${\Gamma }_c$, т.е. от выбора коэффициента расширения вычислительной области $k_c$.

Для численного решения задачи квадратичного программирования (7.8) в настоящей работе использовался алгоритм на основе метода сопряженных градиентов, предложенный в [24] при решении задач одностороннего дискретного контакта для упругой полуплоскости. В [8] этот алгоритм применялся при решении задач одностороннего дискретного контакта для однородной упругой полосы. Следует отметить, что алгоритм [24] позволяет вычислять не только сеточную функцию нормальных напряжений $Q_p\in {ℝ}^L$, но и смещение ${\delta }_2$ и угол поворота ${\phi }_3$ жесткого штампа.

8. Апостериорный анализ численных решений. Для оценки погрешности выполнения сеточных аппроксимаций условий (3.2)$–$(3.5) использовался следующий набор параметров:

${\epsilon }_q=\underset{1\le i\le L}{max}{q}_i/q_{nom}$, ${\epsilon }_f=\left|{\displaystyle \sum _{i=1}^L{\theta }_i{q}_i}-F_2\right|/\left|F_2\right|$, ${\epsilon }_m=\left|{\displaystyle \sum _{i=1}^L{\vartheta }_i{q}_i}-M_3\right|/\left|M_3\right|$ для $M_3\ne 0$

${\epsilon }_{zp}=\underset{1\le i\le L}{min}{z}_i{[q_i\ge 0]}_a/\left|{\delta }_2\right|$, ${\epsilon }_{zm}=\underset{1\le i\le L}{max}{z}_i{[q_i<0]}_a/\left|{\delta }_2\right|$,

где $q_{nom}=F_2/diam{\Gamma }_p$ $–$ средние (номинальные) контактные напряжения; $z_i={\varphi }_i+{\delta }_2+(x_1^i-x_1^c){\phi }_3-w_i$ $–$ конечный зазор (проникание) между упругой полосой и штампом в узле $t_i$ сетки $T_p$; $w_i$ $–$ нормальное перемещение узла $t_i$ сетки $T_p$; ${[\cdot ]}_a$ $–$ скобка Айверсона (функция, равная $1$ для истинного аргумента и равная $0$ в противном случае). Сеточная функция нормальных перемещений $W_p=(w_1,w_2,\dots ,w_L)$, соответствующих решению $Q_p$ задачи квадратичного программирования (7.8), вычислялась с помощью (7.6).

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

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

Также проведено сравнение численных решений рассмотренных ниже задач дискретного контакта для штампов с регулярным поверхностным микрорельефом с решениями, полученными с помощью известного алгоритма [25]. Учитывая возможности последнего, при постановке контактных задач использовалось лишь первое из уравнений равновесия штампа (2.7) и полагалось ${\phi }_3=0$. Выполненные расчеты показали, что при проведении вычислений с двойной точностью среднеквадратичные относительные расхождения решений (сеточных функций нормальных напряжений) не превышают $10^{-7}$.

Получены численные решения ряда задач одностороннего дискретного контакта для стратифицированных упругих полос, в частности, для четырех полос, состоящих из $N=10$ одинаковых слоев, модуль Юнга $E$ и коэффициент Пуассона $\nu$ которых изменялись по одному из следующих законов:

$E({\xi }_2)=E_1+(E_2-E_1){\xi }_2^2$, $\nu ({\xi }_2)={\nu }_1+({\nu }_2-{\nu }_1){\xi }_2$

$E({\xi }_2)=E_1+(E_2-E_1){\xi }_2^4$, $\nu ({\xi }_2)={\nu }_1+({\nu }_2-{\nu }_1){\xi }_2^2$

$E({\xi }_2)=E_2+(E_1-E_2){\xi }_2^2$, $\nu ({\xi }_2)={\nu }_2+({\nu }_1-{\nu }_2){\xi }_2$

$E({\xi }_2)=E_2+(E_1-E_2){\xi }_2^4$, $\nu ({\xi }_2)={\nu }_2+({\nu }_1-{\nu }_2){\xi }_2^2$,

где ${\xi }_2=\left\{N{x}_2/h\right\}$ $–$ «быстрая» координата, $\left\{\cdot \right\}$ $–$ дробная часть числа; $E_1=E^{*}$; $E_2=4E^{*}$; ${\nu }_1=0.4$; ${\nu }_2=0.1$; $E^{*}$ $–$ масштабный множитель (размерная величина, МПа). Графики приведенных выше функций $E({\xi }_2)$ и $\nu ({\xi }_2)$ изображены соответственно на рис. 1,а,б. Номера кривых 1$–$4 соответствуют номерам стратифицированных полос.

Рис. 1

Далее приведен анализ результатов решения задач дискретного контакта для гладких штампов с регулярным поверхностным микрорельефом, состоящим из $K$ микровыступов. Номинальная область контакта полагалась равной ${\Gamma }_p=\{0\le x_1\le d,x_2=h\}$. Форма основания штампов задавалась функцией [8]

$\Phi (x_1)={\Phi }_1(x_1)+{\Phi }_2({\xi }_1)/K$, (9.1)

где ${\Phi }_1(x_1)$ $–$ выпуклая функция, определяющая макроформу штампа; ${\Phi }_2({\xi }_1)$ $–$ строго выпуклая функция, характеризующая форму микровыступов; ${\xi }_1=\left\{K{x}_1/d\right\}$ $–$ «быстрая» координата. Для заданной пары функций ${\Phi }_1(x_1)$ и ${\Phi }_2({\xi }_1)$ формула (9.1) определяет однопараметрическое семейство штампов $\Xi (K)$, в качестве параметра которого выступает число микровыступов $K$. Штампы, принадлежащие к одному семейству, имеют одинаковую макроформу, а их микровыступы являются подобными.

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

${\Phi }_1(x_1)={\beta }_{1}d{(2x_1/d-1)}^{m_1}$, ${\Phi }_2({\xi }_1)={\beta }_{2}d{(2{\xi }_1-1)}^{m_2}$,

где ${\beta }_1\ge 0,{\beta }_2>0,m_1>1,m_2>1$ $–$ безразмерные параметры.

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

$F_2=-f\text{}{E}^{*}d$, $M_3=e{F}_{2}d$,

где $f$ $–$ безразмерный параметр внешней нагрузки; $e$ $–$ безразмерный параметр, характеризующий эксцентриситет равнодействующей внешней нагрузки относительно центра приведения $x^c=(0.5d,h)$.

Необходимое количество $L$ граничных элементов на ${\Gamma }_p$ определялось путем сравнения решений, полученных на вложенных сетках при их двукратном последовательном измельчении. При решении задач для семейства штампов $\Xi (K)$ каждый микровыступ аппроксимировался на сетке из 256 граничных элементов. Наибольшее общее количество элементов сетки для штампа с $K=4096$ микровыступами составило $L=2^{20}=1048576$ элементов.

Установлен ряд закономерностей процесса вдавливания штампов семейства $\Xi (K)$ в упругую полосу. Рассмотрим процесс вдавливания в полосу штампа, имеющего $K$ микровыступов, при изменении параметра внешней нагрузки $f$ в интервале $[0,f^{*}]$ и проведем численный анализ зависимости от этого параметра следующих характеристик контактного взаимодействия: относительной фактической площади контакта

$s\equiv {\displaystyle \underset{{\Gamma }_p}{\int }[{\sigma }_{22}>0]{}_{a}d\text{}{\xi }_1}/d=B_K^1(f)$

осадки штампа ${\delta }_2=B_K^2(f)$; угла поворота штампа ${\phi }_3=B_K^3(f)$.

Введем на $[0,f^{*}]$ равномерную сетку $T_f=\left\{f_m=m\Delta f,\Delta f=f^{*}/M_f\right\}$ и сеточные функции $B_K^i=R_f{B}_K^i,i=1,2,3$, где $R_f$ $–$ оператор ограничения на сетку $T_f$. Для вычисления сеточных функций производилось пошаговое нагружение штампа. При проведении расчетов полагалось, что $M_f=128$.

В ходе вычислительных экспериментов наблюдалась сходимость последовательностей сеточных функций $B_K^i;i=1,2,3$, при увеличении количества микровыступов $K$. В качестве примера для семейства штампов $\Xi (K)$, определяемого набором параметров ( ${\beta }_1=3\cdot 10^{-5}$, ${\beta }_2=10^{-4}$, $m_1=m_2=4$ ), в табл. 1 приведены относительные среднеквадратичные расхождения

${\epsilon }_i={‖B_K^i-B_Q^i‖}_2/{‖B_Q^i‖}_2$; $i=1,2,3$,

где $Q=4096$ $–$ максимальное количество микровыступов в серии расчетов.

Таблица 1. Относительные среднеквадратичные расхождения ${\tilde{\epsilon }}_i=10^5{\epsilon }_i$

Значения параметров внешней нагрузки задавались следующими: $f=5\cdot 10^{-5}$, $e=0.05$. Толщина полосы полагалась равной $h=d$. Коэффициент расширения вычислительной области принимался равным $k_c=16$.

Графики зависимостей $s=B_Q^1(f)$, ${\delta }_2=B_Q^2(f)$ и ${\phi }_3=B_Q^3(f)$ для $Q=4096$ приведены соответственно на рис. 2, а$–$в. Номера кривых 1$–$4 соответствуют номерам стратифицированных упругих полос, параметры которых были указаны выше.

Рис. 2

На основании результатов вычислительных экспериментов можно предположить существование для рассмотренного семейства штампов $\Xi (K)$ предельных кривых $s=B^1(f)$, ${\delta }_2=B^2(f)$ и ${\phi }_3=B^3(f)$.

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

About the authors

A. A. Bobylev

Author for correspondence.
Email: abobylov@gmail.com
Russian Federation, Moscow; Moscow

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1

Download (113KB)

Indexing metadata

3. Fig. 2

Download (251KB)

Indexing metadata