ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ, 2023, том 59, № 1, с.115-129
ЧИСЛЕННЫЕ МЕТОДЫ
УДК 519.624
ЧИСЛЕННОЕ ПОСТРОЕНИЕ ТРАНСФОРМАНТЫ ЯДРА
ИНТЕГРАЛЬНОГО ПРЕДСТАВЛЕНИЯ ОПЕРАТОРА
ПУАНКАРЕ-СТЕКЛОВА ДЛЯ УПРУГОЙ ПОЛОСЫ
© 2023 г. А. А. Бобылев
Рассмотрен оператор Пуанкаре-Стеклова для изотропной стратифицированной упругой
полосы, отображающий на части границы нормальные напряжения в нормальные пере-
мещения. Для построения трансформанты ядра интегрального представления этого опе-
ратора предложен новый подход. Получена вариационная формулировка краевой задачи
для трансформант перемещений. Дано определение и доказаны существование и един-
ственность обобщённого решения задачи. Построен итерационный метод решения вари-
ационных уравнений и на основе принципа сжатых отображений получены условия его
сходимости. Аппроксимация вариационных уравнений проводилась методом конечных
элементов. В результате на каждом шаге итерационного метода требуется решить две
независимые системы линейных алгебраических уравнений, для решения которых приме-
няется метод прогонки. Предложен эвристический алгоритм выбора последовательности
параметров итерационного метода, обеспечивающей его сходимость. Проведена верифика-
ция разработанного вычислительного алгоритма и даны рекомендации по использованию
адаптивных конечно-элементных сеток.
DOI: 10.31857/S0374064123010107, EDN: OCYQEZ
Введение. Операторы Пуанкаре-Стеклова, отображающие на части границы рассматри-
ваемой области граничные условия одного типа в граничные условия другого типа, исполь-
зуются для решения различных классов краевых задач математической физики, в частно-
сти, в методах разделения областей [1, гл. III] и композиции [2, § 17]. Применение операторов
Пуанкаре-Стеклова к решению контактных задач теории упругости с односторонними связями
позволяет получить вариационные формулировки этих задач в виде граничных вариационных
неравенств и задач минимизации граничных функционалов, при численном решении которых
требуется дискретизировать лишь часть границы области - зону возможного контакта, что су-
щественно уменьшает размерность получаемых дискретных задач и снижает вычислительные
затраты [3, 4]. В настоящей работе рассматривается оператор Пуанкаре-Стеклова для изо-
тропной стратифицированной упругой полосы, отображающий на части Γq границы полосы
нормальные напряжения в нормальные перемещения.
Одним из основных подходов к решению плоских задач теории упругости для бесконеч-
ной полосы является применение интегрального преобразования Фурье [5, гл. 1]. В резуль-
тате для оператора Пуанкаре-Стеклова Q : qn → un, отображающего на Γq нормальные
напряжения qn в нормальные перемещения un, можно получить следующее интегральное
представление:
un(x1) = g(x1 - ξ1)qn(ξ1)1,
(1)
Γq
ядро которого имеет вид
1
g(s) =
G(α) cos(αs) dα.
(2)
π
0
Далее, следуя [6], косинус-трансформанту Фурье G(α) ядра g(s) интегрального представле-
ния (1) будем называть передаточной функцией.
115
8
116
БОБЫЛЕВ
В случае однородной полосы выражение для функции G(α) может быть получено ана-
литически [7, § 11]. Для стратифицированной полосы эту функцию удаётся построить с по-
мощью численно-аналитической методики при специальной зависимости её упругих свойств
по толщине, в частности, степенной или экспоненциальной зависимости [8, § 14.5]. В случае
произвольного закона изменения упругих свойств по толщине используются приближённые
подходы, основанные как на замене непрерывно-неоднородной полосы многослойной с кусоч-
но-постоянной зависимостью упругих модулей от координаты [9], так и на прямом численном
интегрировании краевых задач для систем дифференциальных уравнений по поперечной ко-
ординате [6; 10, гл. 2, § 2]. Основные трудности реализации таких подходов обусловлены нали-
чием экспоненциальных составляющих у фундаментальных решений соответствующих систем
дифференциальных уравнений, приводящих к неустойчивости численных процедур решения
задач Коши и их дискретных аналогов и к плохой обусловленности систем линейных алгеб-
раических уравнений (СЛАУ), возникающих при удовлетворении граничных условий. Анализ
публикаций (см., например, [11, 12]) показал, что основным способом преодоления указанных
трудностей является использование метода модулирующих функций, предложенного в [13].
Суть этого метода состоит в выделении в явном виде экспоненциальных составляющих при
построении матрицы Грина, в результате чего проблема сводится к отысканию некоторых
модулирующих функций ограниченной вариации.
В настоящей работе для построения трансформанты ядра интегрального представления
оператора Пуанкаре-Стеклова (передаточной функции G(α)) предложен новый подход, ос-
нованный на использовании вариационной формулировки краевой задачи для трансформант
перемещений. Постановка краевой задачи теории упругости, с помощью решения которой
определяется оператор Пуанкаре-Стеклова, дана в п. 1. Краевая задача для трансформант
перемещений, полученная в результате применения преобразования Фурье, приведена в п. 2.
В п. 3 получена её вариационная формулировка, дано определение и доказаны существование
и единственность обобщённого решения. Для решения полученных вариационных уравнений
в п. 4 построен итерационный метод и на основе принципа сжатых отображений получены
условия его сходимости. Аппроксимация вариационных уравнений проводилась в п. 5 методом
конечных элементов. В п. 6 рассмотрен эвристический алгоритм выбора последовательности
параметров итерационного метода, обеспечивающей его сходимость. Результаты верификации
разработанного вычислительного алгоритма и рекомендации по использованию адаптивных
конечно-элементных сеток приведены в п. 7. В заключении указаны некоторые возможные
обобщения предложенного подхода.
1. Постановка краевой задачи. Сформулируем краевую задачу, с помощью которой
вводится исследуемый оператор Пуанкаре-Стеклова. Пусть изотропная стратифицированная
упругая полоса толщины h занимает область Ω = {x = (x1, x2) R2 : |x1| ∞, 0 x2h}.
Границу полосы x2 0 обозначим Γ0, а границу x2 ≡ h - через Γ1. Под uk(x), σkl(x),
k,l = 1, 2, будем понимать соответственно компоненты вектора перемещений и тензора на-
пряжений в точке x ∈ Ω. Предполагается, что упругая полоса находится в условиях плоской
деформации, деформации малы, а массовые силы и напряжения в недеформированном состо-
янии отсутствуют. Параметры Ламе материала полосы являются произвольными ограничен-
ными функциями координаты x2 : λ = λ(x2) и μ = μ(x2). Требования к гладкости этих
функций будут сформулированы ниже. Из физических соображений следует, что существуют
постоянные λ0 > 0 и μ0 > 0 такие, что выполняются неравенства
λ(x2) λ0, μ(x2) μ0,
0 x2h.
(3)
Напряженно-деформированное состояние полосы Ω описывается системой уравнений
σ11,1 + σ12,2 = 0, σ12,1 + σ22,2 = 0,
σ11 = (λ + 2μ)u1,1 + λu2,2, σ12 = μ(u1,2 + u2,1), σ22 = λu1,1 + (λ + 2μ)u2,2.
(4)
По границе Γ0 полоса соединена с недеформируемым основанием. В случае полного сцеп-
ления граничные условия имеют вид
u1 = u2 = 0 на Γ0.
(5)
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023
ЧИСЛЕННОЕ ПОСТРОЕНИЕ ТРАНСФОРМАНТЫ ЯДРА
117
На части границы Γq Γ1 приложена нормальная нагрузка
σ22 = q, σ21 = 0 на Γq.
(6)
Остальная часть границы Γ1 свободна от внешних нагрузок:
σ22 = σ21 = 0 на Γ1 \ Γq.
(7)
Предполагается, что участок границы Γq является конечным, а главный вектор внешних
усилий имеет ограниченную величину:
∫
diam Γq < ∞,
q dΓq
∞.
(8)
<
Γq
Для завершения постановки краевой задачи необходимо задать условия на бесконечности.
Обычно в качестве таковых используются условия, характеризующие определённый порядок
изменения перемещений и напряжений на бесконечности. Как правило, эти условия носят
чисто математический характер. Более естественным является условие конечности потенци-
альной энергии деформации полосы
[λ(u1,1 + u2,2)2 + 2μ(u21,1 + u22,2) + μ(u1,2 + u2,1)2] dΩ < ∞,
(9)
Ω
которое вполне замыкает постановку задачи и определяет поведение решения на бесконечно-
сти [7, § 1].
2. Применение преобразования Фурье. Введём преобразования Фурье от компонент
перемещений и напряжений обычными формулами
ũk(α,x2) =
uk(x1,x2)exp(-iαx1)dx1,
σkl(α,x2) =
σkl(x1,x2)exp(-iαx1)dx1.
(10)
−∞
-∞
Далее будем полагать, что все преобразования Фурье (10) существуют.
Умножим уравнения (4) на exp(-iαx1) и проинтегрируем по x1 от -∞ до ∞. В ре-
зультате получим систему равенств, которые после интегрирования по частям превращаются
в следующую систему соотношений относительно трансформант
ũk,
σkl и их производных,
рассматриваемых как функции переменной x2 :
iασ11 + σ12 = 0, iασ12 + σ22 = 0,
σ11 =(λ + 2μ)ũ1 + λũ2,
σ12 = μ(ũ1 + iαũ2),
σ22 = iαλũ1 + (λ + 2μ)ũ2.
(11)
Исключив из полученной системы трансформанты напряжений σkl, сделав замену
ũ1(α,x2) = v1(α,x2),
ũ2(α,x2) = iv2(α,x2)
(12)
и положив, что
λ(x2) ∈ C1[0, h], μ(x2) ∈ C1[0, h],
(13)
получим систему двух обыкновенных дифференциальных уравнений
-(μv1) + α(μv2) + αλv2 + α2(λ + 2μ)v1 = 0,
-((λ + 2μ)v2) - α(λv1) - αμv1 + α2μv2 = 0.
(14)
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023
118
БОБЫЛЕВ
Применим далее преобразование Фурье к граничным условиям (5)-(7), предполагая суще-
ствование преобразования Фурье заданной на Γq функции q(x1):
q(α) = q(x1) exp(-iαx1) dx1.
(15)
Γq
В результате, с учётом (11)-(13), получим краевые условия для системы уравнений (14):
v1(α,0) = 0, v2(α,0) = 0,
- (λ(h) + 2μ(h))v2(α, h) - αλ(h)v1(α, h) = p(α), μ(h)v1(α, h) - αμ(h)v2(α, h) = 0,
(16)
где
p(α) =22(α) = iq(α).
(17)
Следовательно, справедлива
Теорема 1. Если существуют преобразования Фурье (10) и (15), а также выполняются
условия (13), то трансформанты компонент перемещений решения краевой задачи (3)-(9)
с учётом (12) удовлетворяют краевой задаче (14), (16).
Определение 1. Решение краевой задачи (14), (16), принадлежащее классу функций
[C2(0, h)]2
[C1[0, h]]2, называют классическим решением.
Замечание 1. Из формул (1), (2), (12), (15)-(17) вытекает соотношение для определения
передаточной функции G(α) с помощью решения краевой задачи (14), (16):
G(α) = -v2(α, h)/p(α).
(18)
Замечание 2. Если α = 0, то рассматриваемая краевая задача (14), (16) распадается на
две независимые краевые задачи относительно функций v1(x2) и v2(x2):
(μv1) = 0, v1(0) = 0, v1(h) = 0,
(19)
((λ + 2μ)v2) = 0, v2(0) = 0, v2(h) = -p(0)/(λ(h) + 2μ(h)).
(20)
Очевидно, что краевая задача (19) имеет тривиальное решение v1 0. Вычислительные ал-
горитмы решения краевых задач типа (20) подробно рассмотрены в литературе по численным
методам (см., например, [14, гл. 9]). Поэтому далее будем считать, что α = 0.
3. Вариационная формулировка. Переход к вариационной формулировке задачи поз-
воляет ослабить требования (13) к гладкости исходных данных(x2), μ(x2)} и определить
обобщённое решение задачи. Введём необходимое для этого пространство вектор-функций
Z = Z × Z, Z = {z ∈ W12(0,h) : z(0) = 0},
(21)
где W12(0, h) - пространство Соболева. Пространство Z является гильбертовым относительно
скалярного произведения
(s, z)Z = (s1, z1)W 1
+ (s2, z2)W 1
,
s = (s1,s2), z = (z1,z2) ∈ Z.
(22)
2
(0,h)
2
(0,h)
Введём в рассмотрение вариационную задачу: для фиксированного α ∈ R найти вектор-
функцию (v1, v2) ∈ Z, удовлетворяющую системе вариационных уравнений
[μv1, w1] - α[μv2, w1] + α[λv2, w1] + α2[(λ + 2μ)v1, w1] = 0 для любой w1 ∈ Z,
[(λ+2μ)v2, w2]+α[λv1, w2][μv1, w2]+α2[μv2, w2]+p(α)w2(h) = 0 для любой w2 ∈ Z, (23)
где [· , ·] - стандартное скалярное произведение в L2(0, h);
λ(x2) ∈ L(0, h), μ(x2) ∈ L(0, h).
(24)
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023
ЧИСЛЕННОЕ ПОСТРОЕНИЕ ТРАНСФОРМАНТЫ ЯДРА
119
Теорема 2. Классическое решение (v1, v2) краевой задачи (14), (16) является решением
вариационной задачи (23).
Доказательство. Пусть w1 и w2 - две произвольные функции из Z. Умножим уравнения
(14) соответственно на w1 и w2, а результаты проинтегрируем по отрезку [0, h]:
h
h
h
h
(μv1)w1 dx2 + α (μv2)w1 dx2 + α λv2w1 dx2
+α2
(λ + 2μ)v1w1 dx2 = 0,
0
0
0
0
h
h
h
h
((λ + 2μ)v2)w2 dx2 - α (λv1)w2 dx2 - α μv1w2 dx2
+α2
μv2w2 dx2 = 0.
0
0
0
0
Первые два интеграла в каждом из полученных уравнений вычислим методом интегрирования
по частям:
h
h
h
h
μv1w1 dx2 - α μv2w1 dx2 + α λv2w1 dx2
+α2
(λ + 2μ)v1w1 dx2 +
0
0
0
0
+ [(-μv1 + αμv2)w1]|h0 = 0,
h
h
h
h
(λ + 2μ)v2w2 dx2 + α λv1w2 dx2 - α μv1w2 dx2
+α2
μv2w2 dx2 -
0
0
0
0
[((λ + 2μ)v2 + αλv1)w2]|h0 = 0.
Учитывая далее краевые условия (16) и принадлежность функций w1 и w2 множеству Z,
получаем вариационные уравнения (23). Теорема доказана.
Определение 2. Решение вариационной задачи (23) называют обобщённым решением
краевой задачи (14), (16).
Используя известные приёмы [15, § 2.14], можно показать, что справедлива
Теорема 3. Решение вариационной задачи (23), если оно существует и обладает вторыми
производными (хотя бы обобщёнными), удовлетворяет (почти всюду) уравнениям (14) и
краевым условиям (16).
Введём на Z билинейные и линейную формы
a1(v,w) = [μv,w] + α2[(λ + 2μ)v,w], a2(v,w) = [(λ + 2μ)v,w] + α2[μv,w],
b1(v,w) = α[μv,w] - α[λv,w], b2(v,w) = α[μv,w] - α[λv,w], l2(w) = -p(α)w(h).
(25)
Лемма 1. Для любого фиксированного α ∈ R при выполнении условий (3) и (24) били-
нейные формы a1(·,·) и a2(·,·) являются непрерывными и положительно определёнными
на Z.
Доказательство. Симметричность билинейных форм очевидна, непрерывность несложно
показать, используя условие (24), а коэрцитивность - с учётом неравенств (3).
Лемма 2. Для любого фиксированного α ∈ R при выполнении условий (24) билинейные
формы b1( · , · ) и b2( · , · ) являются непрерывными на Z.
Доказательство. Непрерывность билинейных форм несложно показать, используя усло-
вие (24), неравенство Коши-Буняковского и неравенство rs (r2 + s2)/2 для любых r, s ∈ R.
Лемма 3. Если существует преобразование Фурье (15), то линейная форма l2(·) непре-
рывна на Z.
Доказательство. Утверждение леммы следует из теоремы вложения C. Л. Соболева
(W12(0, h) ↪→ C[0, h]) [16, § 9].
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023
120
БОБЫЛЕВ
С учётом (25) система вариационных уравнений (23) примет вид
a1(v1,w1) - b1(v2,w1) = 0 для любой w1 ∈ Z,
a2(v2,w2) - b2(v1,w2) - l2(w2) = 0 для любой w2 ∈ Z.
(26)
Образуем далее на Z билинейные и линейную формы
a(v, w) = a1(v1, w1) + a2(v2, w2), b(v, w) = b1(v2, w1) + b2(v1, w2), l(w) = l2(w2).
Используя леммы 1-3, несложно показать, что справедливы следующие утверждения.
Лемма 4. Для любого фиксированного α ∈ R при выполнении условий (3) и (24) били-
нейная форма a(·,·) является непрерывной и положительно определённой на Z.
Лемма 5. Для любого фиксированного α ∈ R при выполнении условий (24) билинейная
форма b( · , · ) является непрерывной на Z.
Лемма 6. Если существует преобразование Фурье (15), то линейная форма l(·) непре-
рывна на Z.
Введём на Z билинейную форму
k(v, w) = a(v, w) - b(v, w)
(27)
и рассмотрим вариационную задачу: для фиксированного α ∈ R найти вектор-функцию v ∈
∈ Z, удовлетворяющую вариационному уравнению
k(v, w) - l(w) = 0 для любой функции w ∈ Z.
(28)
Лемма 7. Вариационное уравнение (28) эквивалентно системе уравнений (26).
Доказательство. Путь v = (v1, v2) - решение системы уравнений (26). Сложив эти урав-
нения, получим с учётом (27) вариационное уравнение (28). Обратно, выбрав в (28) w =
= (w1, 0), получим первое уравнение (26), а выбрав w = (0, w2) - второе уравнение (26).
Лемма доказана.
Разрешимость вариационной задачи (28) зависит от свойств билинейной формы k( · , · ).
Покажем, что эта билинейная форма является положительно определённой на Z.
Из лемм 4 и 5 следует
Лемма 8. Для любого фиксированного α ∈ R при выполнении условий (24) билинейная
форма k( · , · ) является непрерывной на Z.
Непосредственно проверяется, что справедлива
Лемма 9. Билинейная форма k( · , · ) является симметричной на Z.
Для доказательства коэрцитивности билинейной формы k( · , · ) на Z введём гильбертово
пространство вектор-функций Y = [L2(0, h)]4, оснащённое скалярным произведением
y,
y)Y = [ŷ1, y1]+ [ŷ2, y2]+ [ŷ3, y3]+ [ŷ4, y4],
y = (ŷ1, ŷ2, ŷ3, ŷ4),
y = (y1, y2, y3, y4) ∈ Y , (29)
и определим оператор j : Z → Y такой, что jz = (z1, z1, z2, z2) для любого z = (z1, z2) ∈ Z.
Лемма 10. Оператор j : Z → Y является линейным инъективным изометрическим
оператором.
Доказательство. Линейность и инъективность оператора j проверяется непосредствен-
но, а из определений скалярных произведений (22) и (29) следует, что
(jz, jz)Y = (z, z)Z .
Лемма 11. Множество значений R(j) оператора j является собственным замкнутым
линейным подпространством в Y .
Доказательство. Линейность множества R(j) очевидным образом следует из линейно-
сти оператора j, а замкнутость - из полноты пространства Z и изометричности оператора j.
Постоянная вектор-функция (0, 1, 0, 1) ∈ Y не принадлежит R(j), следовательно, R(j) яв-
ляется собственным подмножеством Y . Лемма доказана.
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023
ЧИСЛЕННОЕ ПОСТРОЕНИЕ ТРАНСФОРМАНТЫ ЯДРА
121
Зададим на Y билинейную форму
d(y,y) = (
y,
y)Y ,
(30)
где
α2(λ(x2) + 2μ(x2))
0
0
-αλ(x2)
0
μ(x2) αμ(x2)
0
D(x2, α) =
0
αμ(x2) α2μ(x2)
0
−αλ(x2)
0
0
λ(x2) + 2μ(x2)
– квадратная симметричная матрица-функция четвёртого порядка.
Непосредственно проверяется, что справедлива
Лемма 12. Для любых z,
z ∈ Z выполняется равенство
k(z, ž) = d(j
z,
z).
(31)
Матрицу-функцию D(x2, α) при фиксированных значениях x2 и α можно рассматривать
как квадратную симметричную матрицу четвёртого порядка с постоянными коэффициента-
ми. Современные системы компьютерной алгебры позволяют аналитически вычислить её соб-
ственные значения и соответствующие им собственные векторы. В данной работе с помощью
системы SageMath установлено, что справедлива
Лемма 13. При выполнении условий (3) и (24) матрица D(x2, α) для любого α = 0
почти всюду на отрезке [0,h] имеет четыре вещественных собственных значения:
ω1 = 0, ω2 = [(1 + α2)(λ + 2μ) - κ]/2, ω3 = μ(1 + α2), ω4 = [(1 + α2)(λ + 2μ) + κ]/2,
которым соответствуют ортонормированные собственные векторы
f1 = (0,1,1/α,0)1, f2 = (1,0,0,θ - ϑ)2,
f3 = (0,1,-α,0)3, f4 = (1,0,0 + ϑ)4,
где
κ = [(1 - α2)2(λ + 2μ)2 + 4α2λ2]1/2, θ = (λ + 2μ)(1 - α2)/(2αλ), ϑ = κ/(2αλ),
γ1 = (1 + α2)1/2/α, γ2 = (1 + (θ - ϑ)2)1/2, γ3 = (1 + α2)1/2, γ4 = (1 + (θ + ϑ)2)1/2.
Следствие 1. Для любого α = 0 матрица D(x2, α) является неотрицательно опре-
делённой, а её собственные значения ограничены почти всюду на [0, h] и удовлетворяют
условиям
ω4 > ω3 ω2 > ω1 = 0.
(32)
Следствие 2. Билинейная форма d( · , · ) является неотрицательной на Y .
Из равенства (31) и следствия 2 вытекает, что билинейная форма k( · , · ) также является
неотрицательной на Z. Однако поскольку множество R(j) является собственным подпро-
странством в Y (лемма 11), можно доказать более сильное утверждение - билинейная форма
k(·,·) является коэрцитивной на Z.
Введём квадратную матрицу-функцию F (x2, α) четвёртого порядка, столбцами которой
являются векторы-функции fn(x2, α), n = 1, 4. При фиксированных значениях x2 и α мат-
рицу-функцию F (x2, α) можно рассматривать как квадратную матрицу четвёртого порядка
с постоянными коэффициентами. Из известных результатов линейной алгебры [17, гл. 7] сле-
дует, что справедлива
Лемма 14. Матрица F (x2, α) для любого α = 0 почти всюду на отрезке [0, h] является
ортогональной, т.е.
F-1 = FT и
∥F r∥E4 = ∥r∥E4 для любого r ∈ E4,
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023
122
БОБЫЛЕВ
где E4 - евклидово пространство, и имеет место равенство
FT DF = diag(ω1234).
Следствие 3. Для любого r ∈ E4 и p = FT r выполняются равенства
(Dr, r)E4 = (FT DF p, p)E4 = ω2p22 + ω3p23 + ω4p24.
(33)
Теорема 4. Билинейная форма k( · , · ) является коэрцитивной на Z, т.е. существует
постоянная ϱz > 0 такая, что
k(z, z) ϱz∥z∥2Z для любой z ∈ Z.
Доказательство. Пусть Z1 = {z ∈ Z : ∥z∥Z = 1}, тогда можно положить
ϱz = inf k(z,z).
z∈Z1
Из соотношений (30), (31) и (33) следует, что для любого z ∈ Z1 и q = FT jz имеем
k(z, z) = [ω2q2, q2] + [ω3q3, q3] + [ω4q4, q4].
Из этого равенства и (32) следует, что ϱz 0. Покажем, что ϱz > 0. Доказательство проведём
методом от противного. Пусть ϱz = 0. Тогда, поскольку ω2, ω3, ω4 ограничены почти всюду
на [0, h] (следствие 1), для любого ε > 0 существует z = (z1, z2) ∈ Z1 такой, что для
q =
=FT
z выполняются неравенства
[qn, qn] < ε, n = 2,3,4,
(34)
где
q2 = (z1 + (θ - ϑ)z2)2,
(35)
q3 = (z1 - αz2)3,
(36)
q4 = (z1 + (θ + ϑ)z2)4.
(37)
Из равенств (35) и (37) следует, что
z1 = (θ(γ2q2 - γ4q4) + ϑ(γ2q2 + γ4 q4))/(2ϑ),
z2 = -(γ2 q2 - γ4 q4)/(2ϑ).
Возведём обе части каждого из этих равенств в квадрат и проинтегрируем по отрезку [0, h].
С помощью неравенства Коши-Буняковского и (34) несложно получить оценки
[z1, z1] < ĉ1ε,
[z2, z2] < ĉ4ε,
(38)
где ĉ1 > 0 и ĉ4 > 0 - константы, не зависящие от ε,
q и, следовательно, от z.
Из тождества
x2
z2(x2) = z2(0) +
z2(ξ2)2,
0
учитывая, что z2(0) = 0, для x2 [0, h] с помощью неравенства Коши-Буняковского имеем
неравенство
h
|z2(x2)|
|z2(ξ2)|dξ2 h1/2[z2, z2]1/2,
0
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023
ЧИСЛЕННОЕ ПОСТРОЕНИЕ ТРАНСФОРМАНТЫ ЯДРА
123
из которого с учётом (38) следует оценка
[z2, z2] < ĉ3ε,
ĉ3 = h2ĉ4.
(39)
Далее из соотношения (36) имеем
z1 = γ3 q3 + αz2.
Возведём обе части этого равенства в квадрат и проинтегрируем по отрезку [0, h]. С помощью
неравенства Коши-Буняковского и неравенств (34) и (39) получим оценку
[z1, z1] < ĉ2ε,
(40)
где ĉ2 > 0 - константа, не зависящая от ε,
q и, следовательно, от z.
В результате из (38)-(40) следует, что
∥ z∥2Z < ĉε,
(41)
где
ĉ = ĉ1 + ĉ2 + ĉ3 + ĉ4 > 0 - константа, не зависящая от ε и z. Выбирая в (41) ε < 1/ĉ,
получим неравенство ∥z∥Z < 1, что противоречит условию
z ∈ Z1. Теорема доказана.
Поскольку билинейная форма k( · , · ) является непрерывной и положительно определённой
на Z, а линейная форма l(·) - непрерывной на Z, то из известных результатов (см. [2, гл. 2])
следуют
Теорема 5. Решение вариационной задачи (28) существует и единственно.
Теорема 6. Вариационная задача (28) эквивалентна задаче минимизации функционала
энергии: для фиксированного α ∈ R найти вектор-функцию v ∈ Z такую, что
J (v) = inf {J(w) = k(w, w) - 2l(w)}.
(42)
w∈Z
4. Метод последовательных приближений. Для построения вычислительного алго-
ритма решения рассматриваемой задачи можно использовать метод Бубнова-Галёркина, при-
меняя его непосредственно к решению вариационного уравнения (28), или метод Ритца для
решения задачи минимизации функционала энергии (42). В данной работе используется аль-
тернативный подход - метод последовательных приближений на основе принципа сжатых
отображений.
Лемма 15. Для фиксированного числа τ > 0 существует оператор S : z → s, отобра-
жающий Z в себя, такой, что
a(s, w) = a(z, w) - τ(k(z, w) - l(w)) для любого w ∈ Z.
(43)
Доказательство. Если зафиксировать z, то из лемм 4, 6 и 8 следует, что правая часть
(43) является линейным непрерывным функционалом на Z. Утверждение леммы следует из
теоремы Вишика-Лакса-Мильграма (см. [2, гл. 2, § 14]).
Теорема 7. Вариационная задача (28) эквивалентна задаче нахождения неподвижной
точки v оператора S :
v=Sv.
Доказательство. Пусть v - решение вариационной задачи (28). Положив в (43) z = v,
получим, что s = v. Обратно, пусть v - неподвижная точка v оператора S. Положив в (43)
s = v и z = v, получим, что v удовлетворяет (28). Теорема доказана.
Следствие 4. Неподвижная точка оператора S не зависит от параметра τ > 0.
Получим условия, которым должен удовлетворять параметр τ, чтобы оператор S был
оператором сжатия, т.е. для любых z,
z ∈ Z имело бы место неравенство
∥S z - S
z∥ δ∥z - ž∥,
0 < δ < 1.
(44)
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023
124
БОБЫЛЕВ
Пусть ŝ = S z и š = S ž. Положив в (43) z = z и w = š - ŝ, а затем z = ž и w = ŝ - š,
получим
a(ŝ, š - ŝ) = a(z, š - ŝ) - τ(k(z, š - ŝ) - l(š - ŝ)),
a(š, ŝ - š) = a(ž, ŝ - š) - τ(k(ž, ŝ - š) - l(ŝ - š)).
Сложив эти два равенства, имеем
a(s, s) = a(z, s) - τk(z, s),
(45)
где s = ŝ- š,
z = z - ž.
Билинейная форма a( · , · ) является непрерывной и положительно определённой на Z
(лемма 4), поэтому её можно использовать в качестве скалярного произведения ( · , · )a на Z.
Пространство с новым скалярным произведением обозначим Za. Норма в этом пространстве
будет определяться равенством ∥z∥a = a(z, z)1/2.
Поскольку билинейная форма k( · , · ) является непрерывной и положительно определённой
на Z (леммы 8, 9 и теорема 4), несложно доказать, что справедлива
Лемма 16. Билинейная форма k( · , · ) является непрерывной и положительно определён-
ной на Za, т.е. существуют такие постоянные ϱa > 0 и ρa > 0, что
k(z, z) ϱa∥z∥2a, k(s, z) ρa∥s∥a∥z∥a для любых s, z ∈ Za.
(46)
Из теоремы Вишика-Лакса-Мильграма следует, что существует определённый единствен-
ным образом линейный непрерывный оператор K такой, что
k(s, z) = (Ks, z)a,
∥K∥a = sup
∥Kz∥a ρa.
(47)
∥z∥a=1
Тогда из (45) вытекает, что
∥s∥2a = (z, s)a - τ(K
z, s)a = ((Ia - τK)z, s)a,
(48)
где Ia - единичный (тождественный) оператор в Za.
С помощью неравенства Коши-Буняковского из (48) имеем
∥s∥2a(Ia - τK)z∥a∥s∥a.
(49)
Используя далее (46) и (47), получим оценку
(Ia - τK)z∥2a = ∥z∥2a - 2τ(K
z,z)a + τ2∥Kz∥2a (1 - 2τϱa + τ2ρ2a)∥z∥2a.
(50)
В результате из неравенств (49) и (50) следует, что
∥S
z-Sž∥(1-2τϱa +τ2ρ2a)1/2∥z-
z∥.
(51)
Очевидно, что если τ ∈ (0, 2ϱa2a), то (1 - 2τϱa + τ2ρ2a)1/2 < 1 и, следовательно, выполняется
условие (44). Таким образом, доказана
Теорема 8. Оператор S, определённый формулой (43), является сжимающим для τ из
интервала (0,2ϱa2a).
Для нахождения неподвижной точки v оператора S используется принцип сжатых отоб-
ражений - неподвижная точка отыскивается как предел функциональной последовательности
{vn}, n = 0, 1, . . . , возникающей при решении последовательности вариационных задач
a(vn+1, w) = a(vn, w) - τn+1(k(vn, w) - l(w)) для любого w ∈ Z,
(52)
где τn+1 > 0 - итерационный параметр, который может изменяться от шага к шагу. Начальное
приближение v0 ∈ Z выбирается произвольно.
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023
ЧИСЛЕННОЕ ПОСТРОЕНИЕ ТРАНСФОРМАНТЫ ЯДРА
125
Теорема 9. Если выполняются условия
τ = inf
τn > 0,
τ = supτn < 2ϱa2a,
(53)
n>0
n>0
то итерационный метод (52) сходится при любом начальном приближении v0 ∈ Z к реше-
нию вариационной задачи (28).
Доказательство. Поскольку билинейная форма a( · , · ) является непрерывной и положи-
тельно определённой на Z (лемма 4), а из лемм 4, 6 и 8 вытекает, что правая часть (52) при
фиксированном vn является непрерывной линейной формой, то из теоремы Вишика-Лакса-
Мильграма следует, что решение vn+1 вариационной задачи (52) существует и единственно.
Обозначим через Sn оператор, определённый формулой (43) при τ = τn. Тогда vn+1 = Sn+1vn
и из (51) следуют неравенства
∥vn+1 - vn δ∥vn - vn-1 для n > 0,
(54)
где
δ = max{(1 - 2τϱa + τ2ρ2a)1/2,(1 - 2τϱa + τ2ρ2a)1/2} < 1.
Существование и единственность предела последовательности
{vn} устанавливается с по-
мощью (54) путём рассуждений, аналогичных используемым при доказательстве принципа
сжатых отображений (см. [2, гл. 1, § 4]). Переходя далее в (52) к пределу при n → ∞ и при-
нимая во внимание непрерывность билинейных форм a( · , · ) и k( · , · ), а также условия (53),
получаем утверждение теоремы.
Лемма 17. Вариационное уравнение (52) эквивалентно системе уравнений
a1(vn+11,w1) = (1 - τn+1)a1(vn1,w1) + τn+1b1(vn2,w1) для любого w1 ∈ Z,
a2(vn+12,w2) = (1 - τn+1)a2(vn2,w2) + τn+1b2(vn1,w2) + τn+1l2(w2) для любого w2 ∈ Z.
(55)
Доказательство аналогично доказательству леммы 7.
Система уравнений (55) распадается на два независимых вариационных уравнения, каждое
из которых однозначно разрешимо вследствие непрерывности и положительной определённо-
сти билинейных форм a1( · , · ) и a2( · , · ) на Z. Поэтому несложно показать, что справедлива
Лемма 18. Итерационный метод (52) эквивалентен итерационному методу
a1(vn+1/21,w1) = b1(vn2,w1), a2(vn+1/22,w2) = b2(vn1,w2) + l2(w2) для любых w1,w2 ∈ Z,
(56)
vn+11 = (1 - τn+1)vn1 + τn+1vn+1/21, vn+12 = (1 - τn+1)vn2 + τn+1vn+1/22.
(57)
5. Аппроксимация вариационных задач. Для аппроксимации рассмотренных выше
вариационных задач применим метод конечных элементов. Учитывая определение (21) прост-
ранства Z, для его аппроксимации используем конечномерное подпространство Z ⊂ Z ку-
сочно-линейных функций. Разобъём отрезок [0, h] на M - 1 двухузловых конечных элемента
первого порядка и применим стандартную процедуру метода Бубнова-Галёркина [14, гл. 9,
§ 10] к каждому из вариационных уравнений (56). В результате получим конечномерную ап-
проксимацию итерационного метода (56), (57):
A1Vn+1/21 = B1Vn2, A2Vn+1/22 = B2Vn1 + L2,
(58)
V n+11 = (1 - τn+1)V n1 + τn+1V n+1/21, V n+12 = (1 - τn+1)V n2 + τn+1V n+1/22,
(59)
где A1, A2, B1, B2 - квадратные матрицы порядка M, полученные в результате аппрок-
симации билинейных форм a1( · , · ), a2( · , · ), b1( · , · ) и b2( · , · ) соответственно; L2 RM -
вектор, полученный в результате аппроксимации линейной формы l2(·); Vm1, Vm2 RM , m =
= n,(n + 1), - векторы узловых значений искомых функций на соответствующем шаге итера-
ционного метода; Vn+1/21, Vn+1/22 RM - вспомогательные векторы. Несложно показать, что
B1 = BT2 .
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023
126
БОБЫЛЕВ
Для вычисления элементов матриц A1, A2, B1 и B2 применяются квадратурные форму-
лы. Тип используемых формул определяется конкретным законом изменения упругих свойств
(функций λ(x2) и μ(x2)) по толщине полосы.
На каждом шаге итерационного метода (58), (59) требуется решить две независимые СЛАУ
(58). Поскольку при использовании двухузловых конечных элементов первого порядка матри-
цы A1 и A2 являются трёхдиагональными [14, гл. 9, § 10], то для численного решения СЛАУ
(58) применяется метод прогонки. Несложно показать, что в случае α = 0 обе матрицы A1 и
A2 имеют диагональное преобладание, а если α = 0 - частичное диагональное преобладание.
Следовательно, алгоритм метода прогонки является корректным и устойчивым [18, гл. II].
Применив далее метод Бубнова-Галёркина к вариационному уравнению
(28), полу-
чим СЛАУ
KV = L,
(60)
где
[
]
[
]
[
]
[
]
A1
0
0
B1
V1
0
K =A-B, A=
,
B=
,
V =
,
L=
0
A2
B2
0
V2
L2
Лемма 19. Матрицы A1, A2, A и K являются положительно определёнными.
Доказательство. Поскольку Z ⊂ Z, то положительная определённость матриц A1,
A2, A и K следует из положительной определённости билинейных форм a1(·,·), a2(·,·),
a( · , · ) и k( · , · ) соответственно.
Установим условия, при которых итерационный метод (58), (59) сходится к решению СЛАУ
(60). Непосредственно проверяется, что справедлива
Лемма 20. Итерационный метод (58), (59) эквавалентен итерационному методу
AVn+1 = AVn - τn+1(KVn - L).
(61)
Замечание 3. Итерационный метод (61) представляет собой линейный одношаговый ите-
рационный метод решения СЛАУ (60) с предобуславливателем - матрицей A-1.
Теорема 10 [19, § 34]. Линейный одношаговый итерационный метод (61) сходится к ре-
шению СЛАУ (60), если выполняется условие
sup∥Tn2 < 1,
(62)
n>0
где Tn = I - τnA-1K - матрица перехода n-го шага; I - единичная матрица порядка M;
∥·∥2 - спектральная норма матрицы.
Для определения диапазона значений параметра τn, для которых выполняется условие
(62), используем обобщённое отношение Релея для пучка матриц (K, A):
Φ(V ) = (KV , V )/(AV , V ).
Пусть
ϱ = min
Φ(V ), ρ = max Φ(V ),
(63)
V ∈Q,V =0
V ∈Q,V =0
где Q = Q × Q, Q ⊂ EM - множество узловых значений функций из Z.
Поскольку матрицы A и K являются положительно определёнными, то справедливы
следующие утверждения [20, гл. 15].
Лемма 21. Минимальное значение ϱ обобщённого отношения Релея положительно.
Лемма 22. Спектр матрицы A-1K принадлежит отрезку [ϱ, ρ].
Из теоремы 10 с учётом лемм 21 и 22 следуют известные результаты [19, § 34].
Теорема 11. Линейный одношаговый итерационный метод (61) сходится к решению
СЛАУ (60), если выполняется условие
supτn < 2 = τcr.
n>0
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023
ЧИСЛЕННОЕ ПОСТРОЕНИЕ ТРАНСФОРМАНТЫ ЯДРА
127
Теорема 12. Для линейного стационарного одношагового итерационного метода (61) (ме-
тода простой итерации) асимптотически оптимальным будет значение параметра
τopt = 2/(ϱ + ρ).
(64)
6. Алгоритм выбора последовательности параметровn}. Для вычисления асимп-
тотически оптимального значения τopt в соответствии с формулой (64) требуется найти ми-
нимальное ϱ и максимальное ρ значения обобщённого отношения Релея. Для решения задач
(63) можно использовать степенной метод со сдвигами [20, гл. 15]. Учитывая лемму 21, при
вычислении значения ρ сдвиг полагается равным нулю, а при вычислении ϱ - равным най-
денному значению ρ. Однако вычислительные затраты, оцениваемые по количеству операций
умножения матрицы на вектор, для определения этих двух собственных значений соизмеримы
с трудоемкостью самого итерационного метода (58), (59). Поэтому в данной работе разрабо-
тан эвристический алгоритм выбора последовательности параметровn}, обеспечивающий
сходимость итерационного метода (58), (59).
При вычислении передаточной функции G(α) с помощью формулы (18) на некоторой дис-
кретной сетке значений параметра α в качестве начального приближения для итерационного
метода (58), (59), как правило, используется либо нулевой вектор, либо решение для сосед-
него значения α. Проведённые вычислительные эксперименты по исследованию сходимости
стационарного (τn = τ = const) варианта метода показали, что при таком выборе начального
приближения последовательность разностей vn+12(h) - vn2(h) является знакопостоянной при
τ < τopt - ετ и знакопеременной, если τopt + ετ < τ < τcr, где ετ > 0 - некоторая малая ве-
личина. При τ τcr указанная последовательность также является знакопеременной, однако
возрастающей по абсолютной величине, т.е. итерационный процесс расходится. С учётом этого
предлагается следующий алгоритм выбора последовательности параметровn}.
Пусть определены приближение (vn1, vn2) и значение параметра τn. Используя схему (58),
(59), вычислим два последующих приближения (vn+11, vn+12) и (vn+21, vn+22), полагая τn+2 =
= τn+1 = τn. Далее обозначим
δn+1 = vn+12(h) - vn2(h), δn+2 = vn+22(h) - vn+12(h).
Если выполняется условие
δn+1δn+2 -ωδ2n+1,
(65)
где 0 < ω < 1 - коэффициент сжатия итерационного метода, то два очередных шага счита-
ются удачными, и итерационный процесс продолжается дальше. Если при этом левая часть
неравенства (65) положительна, то значение параметра τ можно увеличить, положив
τn+4 = τn+3 = θτn,
где θ > 1 - коэффициент увеличения параметра τ.
Если условие (65) не выполняется, то необходимо повторить вычисление (n + 1)-го и
(n + 2)-го приближений, положив
τn+2 = τn+1 = ϑτn,
где 0 < ϑ < 1 - коэффициент уменьшения параметра τ.
7. Верификация вычислительного алгоритма. Разработанный вычислительный ал-
горитм построения передаточной функции реализован на языке FORTRAN 2008 в виде пакета
модулей. Для верификации вычислительного алгоритма и программного обеспечения в каче-
стве тестовых выбирались краевые задачи для следующих типов изотропных упругих полос:
- однородные полосы;
- непрерывно-неоднородные полосы, параметры Ламе материала которых являлись экс-
поненциальными функциями λ(x2) = λc exp(κx2) и μ(x2) = μc exp(κx2), где λc, μc и κ -
константы;
- кусочно-однородные полосы с полностью сцепленными слоями.
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023
128
БОБЫЛЕВ
Для упругих полос первого типа известно аналитическое выражение для передаточной
функции [7, § 11], а для остальных строились численно-аналитические решения соответствую-
щих краевых задач [8, § 14]. Анализ этих решений показал, что для рассматриваемого класса
задач характерно следующее поведение решения: с ростом параметра α вблизи точки x2 = h
возникает зона резкого изменения решения, ширина которой уменьшается по мере дальней-
шего увеличения α. При решении таких задач методом конечных элементов целесообразно
использовать адаптивные сетки.
В данной работе применялся следующий алгоритм построения сетки. Задавались общее
количество конечных элементов и размер l(α) наименьшего конечного элемента, одним из
узлов которого является точка x2 = h. Размеры остальных вычислялись таким образом,
чтобы образовывать геометрическую прогрессию по мере удаления элементов от точки x2 =
= h. Эмпирически подобрана зависимость l(α) = 10-3s, где s = 2π/α - “длина волны”,
соответствующая “волновому числу” α.
При наличии разрывов первого рода у функций λ(x2) и μ(x2), а именно такие разрывы
могут быть у модулей упругости реальных слоистых тел, сетка конечных элементов адапти-
ровалась так, чтобы точки разрывов совпадали с узлами конечных элементов. Такой подход
позволяет применять при вычислении элементов матриц A1, A2, B1 и B2 квадратурные
формулы для непрерывных функций.
При использовании адаптивных сеток, состоящих из 200 конечных элементов, среднеквад-
ратичные относительные расхождения численных решений, полученных с помощью разрабо-
танного алгоритма, с указанными выше решениями тестовых задач не превышали 10-5 при
проведении вычислений с двойной точностью.
Для рассмотренного выше алгоритма выбора последовательности параметровn} на ос-
нове проведённых вычислительных экспериментов могут быть рекомендованы для использо-
вания следующие диапазоны значений параметров: ω = 0.4-0.6, θ = 1.3-1.5 и ϑ = 0.4-0.6.
В качестве примера для случая ω = 0.5, θ = 1.4, ϑ = 0.5 в таблице приведена характерная
зависимость общего количества итераций N с учётом повторения итераций при уменьшении
параметра τ от начального значения параметра τ0. Здесь же указаны конечные значения
параметра τN и количество итераций N0 для стационарного итерационного метода (τ = τ0).
В качестве начального приближения решения использовался нулевой вектор. Для данного
расчётного случая также были вычислены τopt = 1.0 и τcr = 1.279. Нетрудно видеть, что
независимо от выбора τ0 в процессе итераций параметр τn стремится к τopt. При отклонении
начального значения τ0 от τopt скорость сходимости нестационарного итерационного метода
снижается, однако не так сильно, как в случае стационарного варианта. При отклонении τ0
от τopt на пять порядков число итераций увеличилось менее чем в четыре раза.
Таблица. Результаты вычислительных экспериментов
τ0
10-5
10-4
10-3
10-2
10-1
100
101
102
103
104
105
τN
0.89
1.21
0.80
1.06
0.72
1.00
0.82
0.72
0.94
1.60
1.43
N
84
68
56
42
30
26
36
44
48
54
62
N0
> 105
> 105
28115
2808
277
26
-
-
-
-
-
В результате вычислительных экспериментов установлено, что при использовании одного
из указанных выше вариантов начального приближения решения и предложенного алгоритма
выбора последовательности параметров
n} итерационный метод (58), (59) сходится для
любого начального значения параметра τ0, изменяется лишь число итераций.
Отметим ещё одно экспериментально установленное свойство итерационного метода (58),
(59). В таблице приведены данные расчётов для 200 конечных элементов. При изменении
числа конечных элементов в диапазоне от 50 до 105 количество итераций как для нестацио-
нарного, так и стационарного вариантов итерационного метода оставалось неизменным. Этот
факт можно объяснить удачным выбором предобуславливателя - матрицы A-1, операция
умножения которой на вектор выполняется путём решения соответствующей СЛАУ методом
прогонки. Кроме того, не изменялись значения τopt и τcr.
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023
ЧИСЛЕННОЕ ПОСТРОЕНИЕ ТРАНСФОРМАНТЫ ЯДРА
129
Заключение. Использование вариационной формулировки краевой задачи для трасфор-
мант перемещений позволило разработать эффективный алгоритм вычисления трансформан-
ты ядра интегрального представления оператора Пуанкаре-Стеклова, отображающего на
части границы изотропной стратифицированной упругой полосы нормальные напряжения в
нормальные перемещения. Разработанный вычислительный алгоритм может быть применён и
в более общем случае - при наличии на границе полосы касательных напряжений. Изменится
лишь вид линейной формы l(w).
Отметим также, что предложенный подход к вычислению передаточной функции может
быть обобщен на случай слоистой полосы при неполном сцеплении слоёв и на случай, когда
упругая стратифицированная полоса сцеплена с упругой однородной полуплоскостью.
СПИСОК ЛИТЕРАТУРЫ
1. Лебедев В.И., Агошков В.И. Операторы Пуанкаре-Стеклова и их приложения в анализе. М., 1983.
2. Лебедев В.И. Функциональный анализ и вычислительная математика. М., 2000.
3. Бобылев А.А. Применение метода сопряжённых градиентов к решению задач дискретного контакта
для упругой полуплоскости // Изв. РАН. Механика твердого тела. 2022. № 2. С. 154-172.
4. Бобылев А.А. Алгоритм решения задач дискретного контакта для упругой полосы // Прикл. ма-
тематика и механика. 2022. Т. 86. № 3. С. 404-423.
5. Уфлянд Я.С. Интегральные преобразования в задачах теории упругости. Л., 1967.
6. Ватульян А.О., Плотников Д.К. К исследованию контактной задачи для неоднородной упругой
полосы // Прикл. математика и механика. 2021. Т. 85. № 3. С. 283-293.
7. Ворович И.И., Александров В.М., Бабешко В.А. Неклассические смешанные задачи теории упруго-
сти. М., 1974.
8. Barber J.R. Contact Mechanics. Cham, 2016.
9. Никишин В.С. Статические контактные задачи для многослойных упругих тел // Механика кон-
тактных взаимодействий. М., 2001. С. 212-233.
10. Айзикович С.М., Александров В.М., Белоконь А.В., Кренев Л.И., Трубчик И.С. Контактные задачи
теории упругости для неоднородных сред. М., 2006.
11. Aizikovich S.M., Galybin A.N., Krenev L.I. Semi-analytical solution for mode I penny-shaped crack in
a soft inhomogeneous layer // Int. J. Sol. Struct. 2015. V. 53. P. 129-137.
12. Trubchik I., Evich L., Ladosha E. Computational model of the deformation of thin gradient coating lying
on nondeformable foundation // AIP Conf. Proc. 2018. V. 1922. Art. 120011.
13. Бабешко В.А., Глушков Е.В., Глушкова Н.В. Методы построения матриц Грина для стратифици-
рованного упругого полупространства // Журн. вычислит. математики и мат. физики. 1987. Т. 27.
№ 1. С. 93-101.
14. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М., 2011.
15. Колтунов М.А., Кравчук А.С., Майборода В.П. Прикладная механика деформируемого твердого
тела. М., 1983.
16. Треногин В.А. Функциональный анализ. М., 2002.
17. Годунов С.К. Современные аспекты линейной алгебры. Новосибирск, 1997.
18. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М., 1978.
19. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М., 1984.
20. Парлетт Б. Симметричная проблема собственных значений. Численные методы. М., 1983.
Московский государственный университет
Поступила в редакцию 18.09.2022 г.
имени М.В. Ломоносова,
После доработки 18.09.2022 г.
Московский центр фундаментальной
Принята к публикации 28.11.2022 г.
и прикладной математики
9
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 59
№1
2023