МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ

имени Н. Э. Баумана
(национальный исследовательский университет)

На правах рукописи

Соколов Андрей Александрович

МАТЕМАТИЧЕСКИЕ МОДЕЛИ НЕЛОКАЛЬНОЙ

ТЕРМОУПРУГОСТИ И ИХ ЧИСЛЕННАЯ
РЕАЛИЗАЦИЯ

1.2.2 — Математическое моделирование, численные методы и комплексы

программ

Диссертация на соискание учёной степени

кандидата физико-математических наук

Научный руководитель:
д.ф.-м.н., доцент
Савельева Инга Юрьевна

Москва — 2024
 2

Оглавление
Стр.

Список сокращений и обозначений . . . . . . . . . . . . . . . . . . . 4

Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

Глава 1. Основные соотношения . . . . . . . . . . . . . . . . . . . . . 16

1.1. Определение нелокального оператора . . . . . . . . . . . . . . . . 16
1.2. Уравнение стационарной теплопроводности . . . . . . . . . . . . . 16
1.3. Уравнение равновесия . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.4. Определение области и функции нелокальности . . . . . . . . . . 19
1.5. Основные результаты и выводы по главе 1 . . . . . . . . . . . . . 25

Глава 2. Численная схема решения на основе метода конечных

элементов . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.1. Общие сведения о методе конечных элементов . . . . . . . . . . . 26
2.2. Аппроксимация уравнений . . . . . . . . . . . . . . . . . . . . . . 27
2.3. Ассемблирование систем уравнений . . . . . . . . . . . . . . . . . 31
2.4. Вычисление производных величин . . . . . . . . . . . . . . . . . . 34
2.5. Основные результаты и выводы по главе 2 . . . . . . . . . . . . . 35

Глава 3. Реализация программного комплекса . . . . . . . . . . . . 36

3.1. Общая структура программного комплекса . . . . . . . . . . . . . 36
3.2. Параллельный алгоритм ассемблирования матриц . . . . . . . . . 42
3.3. Алгоритм аппроксимации области нелокального влияния . . . . . 47
3.4. Оптимизация базисных функций конечных элементов . . . . . . . 47
3.5. Основные результаты и выводы по главе 3 . . . . . . . . . . . . . 51

Глава 4. Анализ результатов расчётов . . . . . . . . . . . . . . . . . 52

4.1. Стратегия исследования и обезразмеривание . . . . . . . . . . . . 52
 3

Стр.

4.2. Основные особенности решений . . . . . . . . . . . . . . . . . . . 53

4.3. Исследование функций нелокального влияния . . . . . . . . . . . 57

4.4. Принципы Сен-Венана и стабильности теплового потока . . . . . 59
4.5. Растяжение пластины со ступенчатым переходом . . . . . . . . . 64
4.6. Задача Кирша с обобщением на эллиптические вырезы . . . . . . 70
4.7. Тепловые деформации в областях с эллиптическими вырезами . . 74
4.8. Основные результаты и выводы по главе 4 . . . . . . . . . . . . . 78

Глава 5. Анализ эффективности программного комплекса

NonLocFEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.1. Тестирование алгоритма ассемблирования матриц . . . . . . . . . 80
5.2. Анализ скорости сходимости при оптимизации базиса конечных
элементов . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.3. Предобуславливание и выбор начального приближения . . . . . . 89
5.4. Основные результаты и выводы по главе 5 . . . . . . . . . . . . . 91

Общие выводы и заключение . . . . . . . . . . . . . . . . . . . . . . . 92

Список литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

Приложение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
 4

Список сокращений и обозначений

Серендиповый конечный элемент: конечный элемент, все узлы кото­

рого находятся на границе элемента.
СЛАУ: Система линейных алгебраических уравнений.
ЭВМ: Электронная вычислительная машина.
MPI: Message Passing Interface.
OpenMP: Open Multi-Processing.
 5

Введение

Задачи термоупругости очень популярны в различных инженерных при­

ложениях, так как температурные деформации могут существенным образом
повлиять на функциональные свойста рассматриваемых объектов, вплоть до их
полного выхода из строя. Особенно популярны такого рода задачи в аэрокосми­
ческой отрасли при моделировании поведения обшивок корпусов и двигателей
летательных аппаратов [33, 92, 118], так как они подвержены очень высоким
и в то же время неравномерным нагружениям [32, 91, 97, 109]. Помимо аэро­
космической отрасли такие задачи могут быть востребованы в строительстве,
особенно в строительстве критической инфраструктуры, такой, например, как
атомные электространции [29, 61]. В гражданской инфраструктуре актуален
анализ влияния таких разрушительных явлений как пожары [59, 100] или обыч­
ная циклическая смена сезона [127, 128], при воздействии которых конструкция
не должна потерять устойчивость. В микроэлектронике эти задачи также не
остаются без внимания [62, 129], а учитывая возрастающее количество вы­
числительных мощностей и популярность микро- и наноэлектромеханических
систем (МЭМС/НЭМС), возникают задачи об эффективном отводе тепловой
энергии [28].
Все перечисленные ранее и многие другие задачи объединяет потребность
в создании новых материалов, которые будут отвечать соответствующим их
использованию требованиям. На сегодняшний день в некоторых отраслях требо­
вания к свойствам материалов становятся уже настолько высокими, что при их
создании необходимо учитывать молекулярную структуру материала на микро-
и наноуровне [35, 126, 132], так как их свойства могут напрямую зависеть
от этого. Такие материалы принято называть структурно-чувствительными, а
создание материалов с наперёд заданными свойствами на сегодняшний день
 6

является одной из сложнейших, но вместе с этим крайне актуальной областью

материаловедения [43].
Вместе с проблемой создания структурно-чувствительных материалов
многие исследователи сталкиваются с проблемой моделирования их поведе­
ния. При рассмотрении наномасштабных структур отсутствует возможность
использовать гипотезу сплошности среды, из-за чего классические модели ме­
ханики сплошной среды не могут даже на качественном уровне передать все
особенности их поведения. Так, например, в наномасштабе перестаёт работать
гипотеза Био — Фурье [76, 137] и наблюдается пониженная чувствительность
к концентраторам напряжений [102, 112]. Этому способствуют такие эффекты,
как микровращения отдельных зёрен материала, микродислокации, различные
дальнодействующие и многие другие масштабные эффекты, которые могут
быть смоделированы только при помощи новых математических моделей.
На сегодняшний день существует большое количество моделей способных
описать различные масштабные эффекты. Однако подходы моделирования мо­
гут достаточно сильно различаться между собой при рассмотрении разных
линейных размеров и временных отрезков. Таким образом, возникают иерархии
моделей, способных качественно и количественно описать различные аспекты
поведения материала на разных масштабах. Это, в свою очередь, приводит к
идее многомасштабного моделирования [17], где, например, некоторые характе­
ристики материала можно вычислять при помощи моделей, находящихся ниже
по иерархии, и передавать полученные в расчётах параметры в вышестоящие
модели или наоборот.
Для механики твёрдого тела одна из возможных иерархий моделей
проиллюстрирована на Рис. 1. Согласно такому представлению модели, исполь­
зующие аппарат квантовой механики [79, 107], находятся на первой ступени
иерархии, их применение ограничено масштабами сопоставимыми с ядрами ато­
мов и простейших молекулярных соединений, состоящих из небольшого числа
 7

атомов, то есть в диапазоне от нескольких ангстрем до нескольких наномет­

ров. На второй ступени иерархии находятся модели молекулярной динамики
[57, 60, 94, 106], такие модели могут описывать поведение сложных соединений,
например, больших полимерных молекул и прочих наномасштабных объектов,
размеры которых не превосходят нескольких десятков нанометров. На третьей
ступени располагаются статистические модели, в частности модели, в осно­
ве которых лежит метод Монте-Карло [111, 136]. В таких моделях расчёты
проводят многократно, а структуру рассматриваемого объекта генерируют слу­
чайным образом по определённым правилам, после чего полученные таким
способом результаты осредняют или вычисляют на их основе вероятностные
характеристики материала. И на последней — четвёртой ступени иерархии сто­
ят континуальные модели, в частности модели механики сплошной среды [34].
Такие модели оперируют гипотезами сплошности среды и абсолютности време­
ни, то есть не учитывают дискретность рассматриваемого вещества.

10 с
Модели
сплошной среды

1с Методы
микромасштаба
1 мс
Молекулярная

1 мкс динамика

1 нс
Квантово-
механическое
моделирование
1 пс
L
1 Å 1 нм 10 нм 100 нм 1 мкм 10 мкм
Рис. 1. Иерархия моделей

Однако у статистических моделей и моделей молекулярной динамики есть

недостаток — анализ объектов при помощи этих моделей без численных экс­
 8

периментов крайне ограничен [98]. Поэтому с середины XX века набирают

популярность модели обобщённой механики сплошной среды, которые распро­
страняют применение моделей высшего уровня на области применения моделей
низшего уровня. Одна из первых таких моделей была предложена в работе
братьев Eugène и François Cosserat [77], где помимо трансляционных степеней
свободы также были учтены и вращательные компоненты движения, которые
связаны с трансляционными рядом соотношений, из-за чего тензор напряже­
ний перестаёт быть симметричным. Позже, спустя пол века, эта теория была
связана с теорией дислокаций в работе V. Günther [96] и дополнена законом
сохранения микроинерции в работе A.C. Eringen [85, 88], в связи с чем тео­
рию начали называть микрополярной теорией упругости. Также к работам, в
которых исследована микрополярная теория упругости, можно отнести рабо­
ты R.D. Mindlin [113, 115, 117], D.B. Bogy [75] и Y.C. Hsu [99]. В них авторы
рассматривали применение этой теории в задачах с концентраторами, возника­
ющим в углах и отверстиях соответсвенно. В то же время теория нашла своё
отражение в работах советских учёных Э.Л. Аэро и Е.В. Кувшинского [21, 22],
а также была рассмотрена в работах Н.Ф. Морозова [45] и Г.Н. Савина [55].
Дальнейшее развитие микрополярной теории упругости привело к появ­
лению микроморфных моделей [44, 86, 135], в которые помимо вращательных
компонент движения могут быть включены дополнительные переменные, свя­
занные с деформацией материала, при этом микрополярная теория упругости
является лишь частным случаем микроморфных моделей. Но стоит учесть тот
факт, что использование таких моделей сопряжено с трудностью определения
материальных коэффициентов.
Список рассмотренных моделей, учитывающих микровращение, а также
авторов, которые занимались их развитием и исследованием, далеко не исчер­
пывающий. Однако стоит также уделить внимание другому классу моделей
обобщённой механики сплошной среды, описывающих дальнодействующие эф­
 9

фекты. Это градиентные модели, которые получили своё развитие в 60-х годах
XX века. Эти модели оперируют высшими производными деформаций, в связи с
чем они и получили такое название. Первые модели градиентной теории упруго­
сти были сформулированы в работах Toupin R.A. [130] и Mindlin R.D. [114, 116],
которые сейчас в литературе принято называть моделями Миндлина — Тупина
[24, 83, 105]. В работе G. Ahmadi и K. Firoozbakhsh [66] эти модели получили
обобщение на температурные деформации. Как и микроморфные, градиентные
модели обладают тем же недостатком — большое количество материальных кон­
стант, которые необходимо определить, поэтому в 90-х годах XX века в работах
E.C. Aifantis и его соавторов [67, 125] была рассмотренна упрощёная модель гра­
диентной теории упругости, в которой напряжения связаны с деформацией и
её вторым градиентом, и, по сравнению с классической моделью, был добавлен
всего один дополнительный материальный параметр внутренней длины.
Существует ещё один класс моделей, описывающих дальнодействующие
эффекты, — это нелокальные модели, которые в отличие от градиентных опе­
рируют интегральными выражениями типа свёртки. Впервые описание таких
моделей было представлено в работе E. Kröner [101], где были рассмотрены
упругие среды с дальнодействующими силами сцепления. Модели нелокаль­
ной упругости в термодинамическом контексте были рассмотрены в работах
D.G.B. Edelen, A.E. Green и N. Laws [81, 82], позже к их работе присоединился
и A.C. Eringen [84, 89]. Исследование условий, обеспечивающих существование
фундаментальных решений было проведено в работе D. Rogula [124]. Вопро­
сы, связанные с существованием и единственностью решений начально-краевых
задач теории нелокальной упругости были рассмотрены в работах S.B. Altan
[69, 71], позже он рассмотрел этот вопрос и для задач нелокальной термоупру­
гости [70, 133]. В конечном итоге, в начале XXI века, A.C. Eringen представил
работу [87], в которой описан единый подход к построению нелокальных теорий
для упругих тел, в связи с чем в литературе нелокальные модели часто назы­
 10

вают моделями Эрингена и их исследованию посвящено достаточно большое

количество работ [64, 90, 131]. При этом между нелокальными и градиентными
моделями существует связь, которая была рассмотрена в работах S.B. Altan и
E.C. Aifantis [68], J. Gao [95] и др.
У нелокальных моделей также есть свои недостатки. Главный из них —
необходимость введения дополнительных условий, так как обычных гранич­
ных условий будет недостаточно [121]. Поэтому в прикладных исследованиях
используют регуляризованные модели, в которых рассматривают комбиниро­
ванные среды, состоящие из локальной и нелокальной фаз. Первые примеры
рассмотрения такого рода сред можно найти в работах C. Polizzoto [122, 123],
где были проанализированы одномерные задачи упругости, а также разрабо­
тан численный метод решения на основе метода конечных элементов. Развитие
этих идей в рамках двумерных задач нелокальной упругости было описано в
работе A.A. Pisano [120].
Анализ моделей нелокальной телопроводности на примере решения одно­
мерных задач был проведён в работе Г.Н. Кувыркина и И.Ю. Савельевой [40].
В это же время Г.Н. Кувыркиным были рассмотрены нелокальные модели тер­
мовязкоупругости в работах [37—39]. Позже в работах И.Ю. Савельевой были
рассмотрены задачи термоудара [51, 54] и вариационные постановки задачи
[50, 52]. Построение термомеханических моделей, оценка тепловых и термо­
упругоих свойств дисперсных структур, а также вариационные формулировки
моделей были представлены в диссертационной работе И.Ю. Савельевой [53]. К
исследованиям нелокальных моделей были подключены и ученики Г.Н. Кувыр­
кина в том числе и автор данной диссертационной работы [41, 42, 103, 104, 110].
Решение задач в нелокальных постановках вызывает достаточно много
сложностей, так как приходится иметь дело с интегро-дифференциальными
уравнениями, которые не всегда имеют аналитические решения даже на про­
стых областях. В этом случае необходимо использовать различные численные
 11

методы, специально адаптированные под данный класс уравнений. В этом

направлении есть уже достаточно большое количество работ, предлагающих
использовать различные методы решения. Наиболее общим и популярным яв­
ляется метод конечных элементов (FEM), который применительно к данному
классу уравнений иногда ещё называют методом нелокальных конечных эле­
ментов (NL-FEM) [120, 123]. Однако его использование сопряжёно с большой
вычислительной сложностью. Поэтому некоторые исследователи используют
его модификацию на основе быстрого преобразования Гаусса (FEMFGT) [74].
Но использование такого подхода сопряжено с проблемами контроля точности.
Для того, чтобы избежать осциляций, необходимо решать задачи на достаточно
подробных сетках, что в некоторых ситуациях лишает данный подход ожидае­
мых преимуществ перед прямым методом решения. Помимо сеточных методов
большой популярностью пользуются и бессеточные подходы на основе радиаль­
ных базисных функций [134], безэлементный метод Галёркина (EFG) и метод
конечных точек (FPM) [65]. Также были предложены подходы с использова­
нием пограничного слоя [64] и на основе полиномов Чебышёва [63], однако на
практике этот метод не применяется в силу своей трудоёмкости, но он может
быть использован для оценки качества решения другими методами.
В рамках текущей работы было принято решение использовать метод
нелокальных конечных элементов, так как данный метод достаточно хорошо
изучен и его легко использовать на областях со сложной геометрией, а большое
количество редакторов и генераторов сеток упрощает процесс моделирования.
Также в работе проведена большая работа по ускорению данного метода, но
чтобы в полной мере реализовать весь потенциал предложенных алгоритмов
был реализован свой собственный програмный комплекс NonLocFEM [56] на
языке программирования C++. Такое решение связано с тем, что многие совре­
менные коммерческие программные комплексы, например Abaqus [1], Ansys [2],
TFlex [16] и др., имеют закрытый программный код. С другой стороны суще­
 12

ствуют открытые программные комплексы, например, [Link] [6], FEniCS [8],

FreeFEM [9] и многие другие. Однако использование открытых программных
комплексов так же может повлечь за собой определённый ряд проблем. Напри­
мер, одной из проблем может стать невозможность эффективно реализовать
тот или иной алгоритм в силу базовых принципов, которые заложены в про­
граммный комплекс. Вторая, наверное даже более серьёзная проблема, связана
с потенциальным конфликтом интересов, так как многие программные комплек­
сы созданы при поддержке зарубежных инвесторов, которые в любой момент
могут ограничить к ним доступ и поэтому важно иметь собственные отечествен­
ные наработки.
Целью исследования является изучение особенностей разработанных дву­
мерных моделей нелокальной теплопроводности и термоупругости, а также
сравнительный анализ решений в случае классических и нелокальных моделей
механики сплошной среды.
Для достижения поставленной цели потребовалось решить задачи.
1. Разработать определяющие соотношения двумерных моделей теплопро­
водности и термоупругости нелокальной среды в интегро-дифференци­
альной форме, а также реализовать эффективные алгоритмы численного
решения на основе метода конечных элементов с последующей реализаци­
ей в виде собственного программного комплекса.
2. Разработать экономичные способы предобуславливания получаемых при
аппроксимации систем линейных алгебраических уравнений (СЛАУ) с це­
лью ускорения сходимости итерационных методов решения.
3. Исследовать особенности нелокальных моделей, сопоставить полученные
результаты в задачах с известными решениями в классической постанов­
ке, определить закономерности.
 13

Научная новизна:
1. Предложены новые эффективные численные алгоритмы для задач нело­
кальной теплопроводности и нелокальной термоупругости на основе мето­
да конечных элементов, которые обладают хорошей масштабируемостью
и предназначены для вычислений на многопроцессорных вычислительных
машинах с общей и распределённой памятью.
2. Разработан собственный программный комплекс NonLocFEM, в котором
реализованы все представленные в работе алгоритмы и методы для моде­
лирования поведения структурно-чувствительных материалов.
3. Получены новые результаты в задачах с известными для классической
постановки решениями, установлены закономерности, свидетельствующие
о снижении роли концентраторов в распределениях полей напряжений и
плотности теплового потока.
4. Исследованы границы спектров собственных чисел матриц и установле­
ны связи между спектрами матриц, ассемблированных в классической и
нелокальной постановках.
Практическая значимость моделей, рассмотренных в диссертации,
состоит в возможности описания поведения термомеханических состояний
структурно-чувствительных материалов. Параметры модели очевидным обра­
зом влияют на решения, что дает возможность точно настраивать модель
для применения на практике. Разработанный программный комплекс, в кото­
ром реализованы численные алгоритмы исследования разработанных моделей,
позволит проводить расчёты на произвольных областях со всеми рассматри­
ваемыми в моделе параметрами, а благодаря открытому исходному коду и
модульной структуре существует возможность редактировать существующие
постановки и добавлять новые типы расчётов при модификации математиче­
ской модели.
 14

Методы исследования. В диссертации использованы как классические

принципы механики деформируемого твёрдого тела, так и новые, относящиеся к
нелокальным теориям теплопроводности и термоупругости, а также численные
методы, в основе которых лежит метод конечных элементов.
Основные положения, выносимые на защиту:
1. Модели нелокальной теплопроводности и термоупругости, позволяющие
описать процессы передачи теплоты и напряжённо-деформированного со­
стояния в структурно-чувствительных материалах.
2. Новые численные алгоритмы решения на основе метода конечных элемен­
тов, адапатированные под многопроцессорные вычислительные системы.
3. Собственный программный комплекс NonLocFEM, в рамках которого ре­
ализованы все рассматриваемые в работе методы решений.
Достоверность результатов гарантирована строгостью и полнотой ис­
пользования возможностей математического аппарата, сравнением результатов
многочисленных проведеннных расчетов с известными аналитическими реше­
ниями и данными, полученными ранее другими авторами.
Апробация работы проводилась в обсуждениях на следующих конфе­
ренциях:
1. Международная научно-техническая конференция «Актуальные пробле­
мы прикладной математики, информатикии и механики» (Воронеж, 2019,
2021);
2. Международная конференция «International Conference of Numerical
Analysis and Applied Mathematics» (Родос, Греция, 2021);
3. Международная научная конференция «Фундаментальные и прикладные
задачи механики» (Москва, 2021);
4. Всероссийская конференция по численным методам решения задач теории
упругости и пластичности (Красноярск, 2023);
 15

5. Международная конференция «Математическое моделирование, числен­

ные методы и инженерное программное обеспечение» (Москва, 2023).
Тема диссертации согласована с тематикой грантов, выделенных на фун­
даментальные исследования.
1. 0705-2020-0047 «Теория дифференциальных уравнений, краевые задачи,
связанные задачи анализа и теории приближений и некоторые их прило­
жения».
2. FSFN-2023-0012 «Разработка математических моделей и методов про­
ектирования изделий ракетно-космической техники из перспективных
конструкционных и функциональных материалов».
3. FSFN-2024-0004 «Разработка математических моделей и методов про­
ектирования изделий ракетно-космической техники из перспективных
конструкционных и функциональных материалов».
Публикации. Основные результаты по теме диссертации изложены
в 5 печатных изданиях, 2 из которых изданы в журналах, рекомендован­
ных ВАК РФ, 3 — в периодических научных журналах, индексируемых Web
of Science и Scopus. Зарегистрирована 1 программа для ЭВМ.
Личный вклад соискателя. Все исследования, представленные в дис­
сертационной работе, а также разработка программного комплекса выполнены
лично соискателем в процессе научной деятельности. Из совместных публи­
каций в диссертацию включен лишь тот материал, который принадлежит
соискателю, заимствованный материал обозначен в работе ссылками.
Объем и структура работы. Диссертация состоит из введения, 5 глав,
заключения и 1 приложения. Полный объём диссертации составляет 111 стра­
ниц, включая 37 рисунков и 9 таблиц. Список литературы содержит 138 на­
именований.
 16

Глава 1. Основные соотношения

1.1. Определение нелокального оператора

Определим линейный интегральный оператор 𝒩 , который представим в

виде взвешенной суммы, где первое слагаемое — это подставляемое в опера­
тор выражение с весовым множителем 𝑝1 , а второе — это же выражение, но
взвешенное по области 𝑆 ′ (𝑥) с весовой функцией φ и весовым параметром 𝑝2 ,
∫︁
𝒩 [𝑓 (𝑥)] = 𝑝1 𝑓 (𝑥) + 𝑝2 φ(𝑥, 𝑥′ )𝑓 (𝑥′ )𝑑𝑆 ′ (𝑥′ ), 𝑥′ ∈ 𝑆 ′ (𝑥). (1.1)
𝑆 ′ (𝑥′ )∩𝑆

Здесь 𝑓 (𝑥) — выражение, описывающее сохраняющуюся физическую субстан­

цию; 𝑝1 > 0 и 𝑝2 ⩾ 0 — весовые параметры модели такие, что 𝑝1 + 𝑝2 = 1;
φ — функция нелокального влияния, нормированная положительная монотон­
но убывающая функция в области 𝑆 ′ (𝑥); 𝑥′ — точка в области 𝑆 ′ (𝑥), в которой
вычисляется влияние на величины находящиеся в точке 𝑥; 𝑆 ′ (𝑥) — область
нелокального влияния с центром в точке 𝑥 ∈ 𝑆; 𝑆 — область занимаемая рас­
сматриваемым телом.
Отметим, что для каждого отдельно взятого физического процесса ℱ
можно определить свой собственный оператор 𝒩ℱ со своим набором весовых
констант 𝑝1 и 𝑝2 , функцией нелокального влияния φ и областью нелокально­
го влияния 𝑆 ′ (𝑥). Однако для упрощения дальнейших выкладок, без потери
общности, ограничимся гипотезой, что для тепловых и механических моделей
параметры нелокальности одинаковые.

1.2. Уравнение стационарной теплопроводности

В произвольной замкнутой области 𝑆 ⊂ R2 с кусочно-гладкой границей

𝜕𝑆 уравнение стационарной теплопроводности имеет вид [34]

∇ · 𝑞 = 𝑞𝑉 , (1.2)
 17

где 𝑞𝑉 — объёмная плотность мощности внутренних источников и стоков тепло­

ты; 𝑞 — вектор плотности теплового потока, который определим как обобщение
гипотезы Био — Фурье [37—39], подставив её в оператор (1.1)
(︁ )︁
𝑞(𝑥) = 𝒩 −λ · ∇𝑇 ,
̂︀ (1.3)

̂︀ = λ𝑖𝑗 𝑒𝑖 ⊗ 𝑒𝑗 — тензор коэффициентов теплопроводности; 𝑇 = 𝑇 (𝑥) —

где λ
поле температуры.
Граничные условия первого, второго и третьего родов для уравнения (1.2)
имеют вид [34]

𝑇 |Γ1 = 𝑇Γ (𝑥), 𝑛 · 𝑞|Γ2 = 𝑓 (𝑥), 𝑛 · 𝑞|Γ3 = α(𝑇𝑎 (𝑥) − 𝑇 (𝑥))), (1.4)

где Γ1 ∪ Γ2 ∪ Γ3 = 𝜕𝑆, Γ1 ∩ Γ2 = Γ1 ∩ Γ3 = Γ2 ∩ Γ3 = ∅; 𝑇Γ (𝑥) и 𝑓 (𝑥) —

функции, задающие температуру и плотность теплового потока на границах Γ1
и Γ2 соответственно; α — коэффициент конвективного теплообмена с внешней
средой; 𝑇𝑎 (𝑥) — температура внешней среды вблизи границы Γ3 . Для простоты
дальнейшего изложения будем предполагать, что функции 𝑇Γ (𝑥), 𝑓 (𝑥) и 𝑇𝑎 (𝑥)
равны нулю во множествах, где они не определены.

1.3. Уравнение равновесия

В произвольной замкнутой области 𝑆 ⊂ R2 с кусочно-гладкой границей

𝜕𝑆 определим уравнение равновесия сплошной среды [34]

∇·σ
̂︀ = 𝑏, (1.5)

̂︀ = σ𝑖𝑗 𝑒𝑖 ⊗ 𝑒𝑗 — тен­
где 𝑏 = 𝑏𝑖 𝑒𝑖 — вектор плотности объёмных сил; σ
зор напряжений. В работе рассматриваем случай несвязанной термоупругой
задачи, поэтому определим тензор напряжений σ
̂︀ через обобщение закона Дю­
амеля — Неймана с использованием оператора (1.1) [37—39]
(︁(︁ )︁)︁
𝑇
σ(𝑥)
̂︀ =𝒩 C
̂︀ · · ̂︀ε−α
̂︀ ∆𝑇 . (1.6)
 18

ε = ε𝑖𝑗 𝑒𝑖 ⊗ 𝑒𝑗 — тензор деформации; C

Здесь ̂︀ ̂︀ = 𝐶𝑖𝑗𝑘𝑙 𝑒𝑖 ⊗ 𝑒𝑗 ⊗ 𝑒𝑘 ⊗ 𝑒𝑙 — тензор

̂︀ 𝑇 = α𝑇𝑖𝑗 𝑒𝑖 ⊗ 𝑒𝑗 — тензор температурных коэф­

коэффициентов упругости; α
фициентов линейного расширения; ∆𝑇 = 𝑇 − 𝑇0 — разница между текущим
распределением температуры 𝑇 и распределением 𝑇0 при котором отсутствуют
температурные деформации.
Далее будем считать, что тело является линейно-упругим и изотропным.
В случае плоского напряжённого состоянии, компоненты тензора упругости C
̂︀

будут определены следующим образом [34]

ν𝐸 𝐸
𝐶𝑖𝑗𝑘𝑙 = 2
δ𝑖𝑗 δ𝑘𝑙 + (δ𝑖𝑘 δ𝑗𝑙 + δ𝑖𝑙 δ𝑗𝑘 ),
1−ν 2(1 + ν)

где 𝐸 — модуль Юнга; ν — коэффициент Пуассона; δ𝑖𝑗 — дельта Кронекера.

Если же рассмотрен случай плоского деформированного состояния, то компо­
ненты тензора упругости имеют аналогичную форму записи

̃︀𝐸
ν ̃︀ 𝐸̃︀
𝐶𝑖𝑗𝑘𝑙 = δ δ
𝑖𝑗 𝑘𝑙 + (δ𝑖𝑘 δ𝑗𝑙 + δ𝑖𝑙 δ𝑗𝑘 ),
1−ν ̃︀2 2(1 + ν)
̃︀

̃︀ = 𝐸/(1 − ν2 ) и ν
однако, здесь 𝐸 ̃︀ = ν/(1 − ν). Также будем считать, что тело
расширяется равнонаправлено, поэтому тензор температурных коэффициентов
линейного расширения будет диагональным и иметь всего один коэффициент
α𝑇 , то есть [34]

̂︀ 𝑇 = α𝑇 ̂︀I2 .
α

Примем гипотезу, что деформации достаточно малы, поэтому для опре­

ε воспользуемся соотношениями
деления компонент тензора деформации ̂︀
Коши [34]

∇𝑢 + (∇𝑢)𝑇 𝑢𝑖,𝑗 + 𝑢𝑗,𝑖

ε=
̂︀ = 𝑒𝑖 ⊗ 𝑒𝑗 ,
2 2

где 𝑢 — вектор перемещения.

 19

Будем рассматривать граничные условия первого и второго родов [34],

также именуемые кинематическими и силовыми соответственно,

𝑢|Γ4 = 𝑑(𝑥), 𝑛 · σ|
̂︀ Γ5 = 𝑝(𝑥), (1.7)

где 𝑑(𝑥) = 𝑑𝑖 (𝑥)𝑒𝑖 — вектор перемещений на границе Γ4 ; 𝑝(𝑥) = 𝑝𝑖 (𝑥)𝑒𝑖 —

вектор плотности поверхностностного нагружения на границе Γ5 . Помимо это­
го будем рассматривать комбинированные граничные условия, когда по одной
компоненте задано перемещение, а по другой поверхностное нагружение. Как и
в случае с граничными условиями уравнения теплопроводности, для простоты
будем считать, что функции задающие граничные условия уравнения равнове­
сия (1.6) будут равны нулю на границах, на которых они не определены.

1.4. Определение области и функции нелокальности

В определении оператора (1.1) нет ограничения на выбор области нело­

кального влияния 𝑆 ′ (𝑥). Она может быть как неограниченной и включать в
себя всю расчётную область, так и замкнутой, покрывая лишь часть рассмат­
риваемого тела. В любом случае, выбор области 𝑆 ′ (𝑥) подразумевает так же и
выбор функции нелокального влияния φ. С практической точки зрения, следу­
ет выбирать такую функцию φ, чтобы интеграл от неё по области 𝑆 ′ (𝑥) был
в рамках заданной точности близким к единице [85]. Вместе с этим область
𝑆 ′ (𝑥) должна быть достаточной для аппроксимации наблюдаемых явлений, но
в то же время, она не должна покрывать всю область занимаемую телом, так
как на практике при аппроксимации уравнений это позволит использовать раз­
реженные матрицы для хранения коэффициентов СЛАУ [47] и значительно
облегчит численные расчёты, повысив общую эффективность использования
вычислительных ресурсов.
Выбор области нелокального влияния 𝑆 ′ (𝑥) является нетривиальной за­
дачей, где в первую очередь стоит опираться на структуру рассматриваемого
 20

материала [85]. Поэтому рассмотрим наиболее общий (пусть и не исчерпываю­

щий) случай и представим 𝑆 ′ (𝑥) в виде фигуры ограниченной кривой Ламэ [49],
изображённой на Рис. 1.1 при различных параметрах 𝑛 > 0. У такого семейства
фигур есть также параметры, отвечающие за длины главных полуосей 𝑟1 > 0
и 𝑟2 > 0. На основе этих параметров можем определить метрическую функцию

(︂⃒ ′
⃒𝑛 ⃒ ′
⃒𝑛 )︂ 1
⃒ 𝑥1 − 𝑥1 ⃒
ρ𝑛 (𝑥, 𝑥′ ) = ⃒⃒ ⃒ + ⃒ 𝑥2 − 𝑥2 ⃒ 𝑛 ,
⃒ ⃒
(1.8)
𝑟1 ⃒ ⃒ 𝑟2 ⃒

использовав которую приступим к построению всевозможных функций нело­

кального влияния φ, соблюдая правило, что функция φ должна монотонно
убывать по мере роста функции ρ𝑛 .

x2
1.0

0.5 n = 0.5
n=1
x1 n=2
-2 -1 1 2
n=5
-0.5

-1.0
Рис. 1.1. Кривые Ламэ при различных параметрах 𝑛

Воспользовавшись метрической функцией ρ𝑛 (1.8), построим семейство

полиномиальных функций нелокального влияния, определённых в ограничен­
ной области [93],
⎧
⎨𝐴(1 − ρ𝑛 (𝑥, 𝑥′ )𝑝 )𝑞 , ρ𝑛 (𝑥, 𝑥′ ) ⩽ 1,
⎪
φ𝑃𝑝,𝑞 (𝑥, 𝑥′ ) = (1.9)
⎩0, ′
ρ𝑛 (𝑥, 𝑥 ) > 1,
⎪

где 𝑝 > 0 и 𝑞 > 0 — параметры управляющие плотностью распределения

функции; 𝐴 — нормирующий множитель. Для определения нормирующего мно­
жителя 𝐴 проведём нормировку области 𝑆 ′ (𝑥) вдоль каждой из осей таким
 21

образом, чтобы исключить параметры 𝑟1 и 𝑟2 из метрической функции ρ𝑛 (1.8).

𝑥′ ) введём аналог полярной
Далее на получившейся безразмерной области 𝑆̃︀′ (̃︀
системы координат с обобщением на параметр 𝑛, где координату радиуса ρ и
угла θ вычислим по следующим формулам
⎧
1
𝑥1 |𝑛 + |̃︀
𝑥2 |𝑛 ) 𝑛 ,
⎪
⎨ρ = (|̃︀
⎪
(︂ )︂
𝑥
̃︀2
⎩θ = arctg .
⎪
⎪
𝑥
̃︀1
Тогда обратная зависимость координат принимает вид
⎧
ρ
⎨𝑥̃︀1 = 1 ,
⎪
⎪
⎪
(1 + tg𝑛 (θ)) 𝑛
ρ tg(θ)
𝑥 = 1 ,
⎪
⎪
⎪ ̃︀2
(1 + tg𝑛 (θ)) 𝑛
⎩

на основе которой можем вычислить якобиан необходимый для интегрирования

внутри области 𝑆̃︀′ (𝑥)
⎛ ⎞
𝜕̃︀
𝑥1 𝜕̃︀
𝑥1
⎜ 𝜕ρ 𝜕θ ⎟ ρ
𝐽𝑛 = det ⎜
⎝ 𝜕̃︀
𝑥2 𝜕̃︀
⎟=
𝑥2 ⎠ cos2 θ (1 + tg𝑛 (θ)) 𝑛2 .
𝜕ρ 𝜕θ
Теперь проинтегрируем по области 𝑆̃︀′ (𝑥) полиномиальную функцию нелокаль­
ного влияния φ𝑃𝑝,𝑞 (1.9)
∫︁2π ∫︁1
𝐴(1 − ρ𝑝𝑛 )𝑞 𝐽𝑛 𝑑ρ𝑑θ = 1,
0 0
откуда, перейдя обратно к размерной области 𝑆 ′ (𝑥), можем установить, что
величина нормировочного параметра 𝐴 равна выражению
𝑛𝑝
𝐴= (︂)︂ (︂ )︂ ,
1 1 2
4𝑟1 𝑟2 B , B ,𝑞 + 1
𝑛 𝑛 𝑝
где B — бета функция Эйлера [18]. Отдельно отметим, что при стремлении
параметра 𝑛 к бесконечности, нормировочный множитель принимает значение
𝑝
𝐴= (︂ )︂ ,
2
8𝑟2 𝑟2 B ,𝑞 + 1
𝑝
 22

а метрическая функция ρ𝑛 (1.8) вырождается в следующую

′⃒ ⃒ ′⃒
(︂⃒ ⃒ ⃒ ⃒)︂
𝑥 1 − 𝑥 1⃒ ⃒ 2𝑥 − 𝑥
ρ∞ (𝑥, 𝑥′ ) = max ⃒⃒ 2⃒
⃒
,⃒ . (1.10)
𝑟1 ⃒ 𝑟2 ⃒
Далее при использовании данного семейства функций, в случае когда длины
полуосей 𝑟1 и 𝑟2 равны, будем обозначать их одним символом 𝑟 и называть его
радиусом нелокальности.
На Рис. 1.2 представлены распределения полиномиальных функций нело­
кального влияния (1.9) в сечении вдоль оси 𝑥1 , где показано, что увеличение
параметра 𝑝 делает распределение функции более равномерным и в пределе та­
кое распределение стремится к константе обратно пропорциональной площади
заключённой в область 𝑆 ′ (𝑥). Увеличение параметра 𝑞 концентрирует распреде­
ление в центре области и в пределе распределение стремится к дельта-функции
Дирака.
P P
φp,q φp,q
1.0
q=1 p=1
p=1 3.0
r1 = 1 0.8 r1 = 1
2.5
x=0 p=2 x=0 q=3
0.6 2.0

p=5 1.5 q=2

0.4
1.0
0.2 p = 50
0.5 q=1
x1′ x1′
-1.0 -0.5 0.5 1.0 -1.0 -0.5 0.5 1.0

а) б)
Рис. 1.2. Портреты полиномиальных функций влияния в сечении вдоль оси 𝑥′1
при вариации (а) параметра 𝑝 и (б) параметра 𝑞

Аналогично можем определить семейство экспоненциальных функций,

для которых область влияния бесконечная

′ ′ 𝑝
φ𝐸
𝑝,𝑞 (𝑥, 𝑥 ) = 𝐴 exp (−𝑞ρ𝑛 (𝑥, 𝑥 ) ) . (1.11)

Параметры 𝑝 > 0 и 𝑞 > 0 — параметры плотности распределения; 𝐴 — нормиро­

вочный коэффициент. Для определения параметра 𝐴 проделаем аналогичную
 23

процедуру, которая была рассмотрена в полиномиальном семействе функций,

но область интегрирования неограничена, поэтому при переходе в новую систе­
му координат необходимо вычислить несобственный интеграл
∫︁2π ∫︁∞
𝐴 exp (−𝑞ρ𝑝 ) 𝐽𝑛 𝑑ρ𝑑θ = 1,
0 0

откуда находим, что значение нормировочного множителя 𝐴 равно следующе­

му выражению
1 2
4 𝑛 𝑛𝑝𝑞 𝑝
𝐴= (︂ )︂ (︂ )︂ ,
1 1 2
8𝑟1 𝑟2 B , Γ
2 𝑛 𝑝

а при стремлении 𝑛 к бесконечности он принимает следующую форму

2
𝑝𝑞 𝑝
𝐴= (︂ )︂ ,
2
8𝑟1 𝑟2 Γ
𝑝

где Γ — гамма-функция [18]. Обратим внимание, что при 𝑞 = 0.5, 𝑝 = 2 и

𝑛 = 2 получаем функцию нормального распределения Гаусса, для которой па­
раметры 𝑟1 и 𝑟2 можем определить по правилу «3 сигма» [46]. Здесь, как и в
случае с полиномиальным семейством функций, в случае равенства 𝑟1 и 𝑟2 , бу­
дем обозначать их одним символом 𝑟, но называть его будем дисперсионным
параметром нелокальности.
На Рис. 1.3 представлены распределения экспоненциальных функций
нелокального влияния (1.11) в сечениях вдоль оси 𝑥1 . Здесь параметры 𝑝 и
𝑞 имеют тот же смысл, что и у полиномиального семейства.
Для экспоненциального семейства функций при заданой ограниченной об­
ласти 𝑆 ′ (𝑥) возникает потребность в определении параметров 𝑟1 и 𝑟2 таких,
чтобы нормировка функции была в рамках заданного квантиля 0 < 𝑄 < 1.
Для этого снова применим процедуру перехода к обобщённой полярной систе­
ме координат, исключим параметры 𝑟1 и 𝑟2 из метрической функции (1.8) путём
 24

E E
φp,q φp,q
q=1 0.35 p = 10 p=1 1.4
r1 = 1 0.30 r1 = 1 1.2
x=0 0.25 p=2 x=0 1.0 q = 3
0.20 0.8
p = 1.3
0.15 0.6
q=2
0.10 p=1 0.4
0.05 0.2 q=1
x1′ x1′
-3 -2 -1 1 2 3 -3 -2 -1 1 2 3

а) б)
Рис. 1.3. Портреты экспоненциальных функций влияния в сечении вдоль оси
𝑥′1 при вариации (а) параметра 𝑝 и (б) параметра 𝑞

её обезразмеривания, тогда верхний предел интегрирования будем считать пе­

ременным и равным 𝑅

∫︁2π ∫︁𝑅
𝐴 exp (−𝑞ρ𝑝 ) 𝐽𝑛 𝑑ρ𝑑θ = 𝑄.
0 0

После интегрирования приходим к уравнению

(︂ )︂
2
Γ , 𝑞𝑅𝑝
𝑝
1− (︂ )︂ = 𝑄, (1.12)
2
Γ
𝑝
из которого путём численного решения можем найти длину 𝑅 при заданных
параметрах 𝑝, 𝑞 и квантиля 𝑄. Отметим, что данное уравнение не зависит от
параметра 𝑛, что упрощает анализ при выборе функций влияния. Кроме то­
го в практических расчётах значение параметра 𝑄 следует выбирать достато
близким к 1. Так, например, для функции нормального распределения, при
условии 𝑄 = 0.99 получим значение 𝑅 = 3.03485, что соответствует известному
правилу «3 сигма» [46]. После нахождения параметра 𝑅, значения параметров
𝑟1 и 𝑟2 можем найти просто поделив соответствующую длину полуоси области
𝑆 ′ (𝑥) на величину 𝑅.
 25

Отметим, что параметр 𝑞 в экспоненциальном семействе функций (1.11)

является избыточным, так как его вариация напрямую связана с параметрами
𝑟1 и 𝑟2 . При подборе параметров 𝑟1 и 𝑟2 по вышеизложенному алгоритму при раз­
личных 𝑞, распределения функций будут одинаковыми. Однако этот параметр
введён намеренно, так как он упрощает управление распределением функций и
с некоторыми оговорками будет использован в дальнейших исследованиях при
сравнении с полиномиальным семейством функций.

1.5. Основные результаты и выводы по главе 1

1. Определён интегральный нелокальный оператор; с его помощью опреде­

лены уравнения стационарной теплопроводности и равновесия в нелокаль­
ных постановках, которые представлены в интегро-дифференциальной
форме.
2. Предложены два семейства функций нелокального влияния: полиноми­
альное семейство функций, с ограниченной областью определения, и
экспоненциальное семейство функций, с неограниченной областью опре­
деления.
 26

Глава 2. Численная схема решения на основе метода конечных

элементов

2.1. Общие сведения о методе конечных элементов

В качестве численного метода решения для уравнений (1.2) и (1.5)

выберем метод конечных элементов с использованием изопараметрических
конечных элементов [72, 138]. Для этого на области 𝑆 введём сетку конечно­
элементной модели 𝑆ℎ , которая включает в себя множества номеров узлов и
связей между ними, образующих непосредственно сами элементы. Каждый эле­
мент 𝑒 ∈ 𝑆ℎ содержит в себе множества узлов {𝑥𝑖 }𝑖∈𝐼 𝑒 и базисных функций
{𝑁𝑖𝑒 }𝑖∈𝐼 𝑒 таких, что

𝑁𝑖𝑒 (𝑥𝑗 ) = δ𝑖𝑗 , 𝑖,𝑗 ∈ 𝐼 𝑒 ,

𝑁𝑖𝑒 (𝑥) = 0, 𝑥 ∈
/ 𝑆 𝑒,
∑︁
𝑁𝑖𝑒 (𝑥) = 1, 𝑥 ∈ 𝑆 𝑒 ,
𝑖∈𝐼 𝑒

где 𝐼 𝑒 — множество индексов узлов элемента 𝑒; δ𝑖𝑗 — дельта Кронекера; 𝑥𝑗 — зна­

чения глобальных координат в узлах сетки; 𝑆 𝑒 — область элемента 𝑒.
Для каждого конечного элемента 𝑒 введём локальную систему коорди­
нат 𝑂ξ𝑒1 ξ𝑒2 . Отображение из локальной системы координат 𝑂ξ𝑒1 ξ𝑒2 в глобальную
𝑂x1 x2 будем строить следующим образом

𝑥(ξ𝑒 ) = 𝑁𝑖𝑒 (ξ𝑒 ) 𝑥𝑖 , 𝑖 ∈ 𝐼 𝑒 , 𝑒 ∈ 𝑆ℎ .

Тогда матрицу Якоби перехода из локальной системы координат в глобальную

можем представить в виде
(︂ 𝑒 )︂ (︂ )︂−1 (︂ 𝑒 −1
)︂
𝑒 𝜕ξ 𝜕𝑥 𝜕𝑁
J
̂︀ = = 𝑒 ≈ 𝑥𝑖 𝑒𝑖 ,
𝜕𝑥 𝜕ξ 𝜕ξ
 27

а вычисление производных функций форм относительно глобальных координат

примет следующий вид

𝜕𝑁𝑖𝑒 𝜕𝑁𝑖𝑒 𝜕ξ𝑒𝑗

= .
𝜕𝑥𝑘 𝜕ξ𝑒𝑗 𝜕𝑥𝑘

Аппроксимированную границу Γℎ ⊂ 𝑆ℎ представим в виде набора од­

номерных элементов, располагающихся на гранях двумерных элементов. При
интегрировании внешних воздействий на границах области также возникает
необходимость в аппроксимации якобиана, формулу которого можно опреде­
лить в следующем виде
⎯
⎸ 2 (︂ 𝑒 2
)︂
𝑒
⎸∑︁ 𝜕𝑁 𝑖
𝐽 =⎷ 𝑥𝑗𝑖 𝑒 , 𝑖 ∈ 𝐼 𝑒 , 𝑒 ∈ Γℎ .
𝑗=1
𝜕ξ

2.2. Аппроксимация уравнений

Спроецируем уравнения (1.2) и (1.5) на функцию 𝑁𝑛𝑒 , где 𝑛 ∈ 𝐼 𝑒 , 𝑒 ∈ 𝑆ℎ

∫︁
𝑁𝑛𝑒 (∇ · 𝑞 − 𝑞𝑉 ) 𝑑𝑆 = 0,
𝑆
∫︁
𝑁𝑛𝑒 (∇ · σ
̂︀ − 𝑏)𝑑𝑆 = 0.
𝑆

Проинтегрируем по частям первое слагаемое каждого уравнения, тогда пользу­

ясь формулой Грина получаем следующие равенства
∫︁ ∫︁ ∫︁
∇𝑁𝑛 · 𝑞𝑑𝑆 − 𝑛 · 𝑞𝑑Γ = 𝑁𝑛𝑒 𝑞𝑉 𝑑𝑆,
𝑒

𝑆 𝜕𝑆 𝑆
∫︁ ∫︁ ∫︁
∇𝑁𝑛𝑒 · σ𝑑𝑆
̂︀ − 𝑛 · σ𝑑Γ
̂︀ = 𝑁𝑛𝑒 𝑏𝑑𝑆.
𝑆 𝜕𝑆 𝑆

Подставим определения граничных условий в уравнения теплопроводности

(1.4) и равновесия (1.7). Интегралы по границе 𝜕𝑆 разбиваем на суммы ин­
 28

тегралов
∫︁ ∫︁ ∫︁ ∫︁ ∫︁
∇𝑁𝑛 · 𝑞𝑑𝑆 + α𝑁𝑛 𝑇 𝑑Γ = 𝑁𝑛 𝑞𝑉 𝑑𝑆 + 𝑁𝑛 𝑓 𝑑Γ + α𝑁𝑛𝑒 𝑇𝑎 𝑑Γ,
𝑒 𝑒 𝑒 𝑒

𝑆
∫︁Γ3 𝑆
∫︁ Γ2
∫︁ Γ3

∇𝑁𝑛𝑒 · σ𝑑𝑆
̂︀ = 𝑁𝑛𝑒 𝑏𝑑𝑆 + 𝑁𝑛𝑒 𝑝𝑑Γ.
𝑆 𝑆 Γ5

Воспользовавшись определениями вектора плотности теплового потока (1.3) и

тензора напряжений (1.6), а также определением оператора (1.1), приходим к
следующим равенствам относящимся к уравнению теплопроводности
⎛ ⎞
∫︁ ∫︁
∇𝑁𝑛𝑒 · ⎝−𝑝1 λ
̂︀ · ∇𝑇 𝑑𝑆 − 𝑝2 φ(𝑥, 𝑥′ )λ
̂︀ · ∇𝑇 𝑑𝑆 ′ (𝑥)⎟
⎠ 𝑑𝑆+
⎜

𝑆 𝑆 ′ (𝑥)∩𝑆
∫︁ ∫︁ ∫︁ ∫︁
+ α𝑁𝑛𝑒 𝑇 𝑑Γ = 𝑁𝑛𝑒 𝑞𝑉 𝑑𝑆 + 𝑁𝑛𝑒 𝑓 𝑑Γ + α𝑁𝑛𝑒 𝑇𝑎 𝑑Γ, (2.1)
Γ3 𝑆 Γ2 Γ3

и уравнению равновесия
∫︁ (︃
(︁ )︁
𝑒 𝑇
∇𝑁 · 𝑝1 C · · ̂︀
𝑛
̂︀ ε−α
̂︀ ∆𝑇 𝑑𝑆+
𝑆
∫︁ )︃
(︁ )︁
+ 𝑝2 φ(𝑥, 𝑥′ )C
̂︀ · · ̂︀ ̂︀ 𝑇 ∆𝑇 𝑑𝑆 ′ (𝑥) 𝑑𝑆 =
ε−α
𝑆 ′ (𝑥)∩𝑆
∫︁ ∫︁
= 𝑁𝑛𝑒 𝑏𝑑𝑆 + 𝑁𝑛𝑒 𝑝𝑑Γ. (2.2)
𝑆 Γ5

Аппроксимируем температуру 𝑇 и вектор перемещения 𝑢 на элементе

𝑒 𝑒
𝑇 (𝑥) = 𝑇𝑚 𝑁𝑚 (𝑥), 𝑢(𝑥) = 𝑢𝑚 𝑁𝑚 (𝑥), 𝑥 ∈ 𝑆 𝑒,

где 𝑇𝑚 и 𝑢𝑚 — искомые значения температуры и вектора перемещения в узле

𝑚 ∈ 𝐼 𝑒 . Тогда аппроксимация градиента температуры и тензора деформации
ε примут следующий вид
̂︀

𝑒
∇𝑇 = 𝑇𝑚 𝑁𝑚,𝑘 𝑒𝑘 , (2.3)
1 (︀ 𝑒 𝑒 𝑇
)︀ 1 𝑒 𝑒
ε=
̂︀ 𝑢𝑚 ∇𝑁𝑚 + (𝑢𝑚 ∇𝑁𝑚 ) = (𝑢𝑚𝑘 𝑁𝑚,𝑙 + 𝑢𝑚𝑙 𝑁𝑚,𝑘 )𝑒𝑘 ⊗ 𝑒𝑙 . (2.4)
2 2
 29

Подставим аппроксимированные значения (2.3) и (2.4) в уравнения (2.1) и (2.2)

соответственно и перейдём к индексной форме записи. Разделив локальные и
нелокальные слагаемые получим системы уравнений для уравнения теплопро­
водности
∫︁ ∫︁ ∫︁
′
− 𝑝1 𝑇𝑚 𝑒
λ𝑖𝑗 𝑁𝑛,𝑖 𝑒
𝑁𝑚,𝑗 𝑑𝑆 − 𝑝2 𝑇𝑚′ 𝑒
𝑁𝑛,𝑖 φ(𝑥, 𝑥′ )λ𝑖𝑗 𝑁𝑚
𝑒 ′
′ ,𝑗 𝑑𝑆 (𝑥)𝑑𝑆+

𝑆 𝑆 𝑆 ′ (𝑥′ )∩𝑆
∫︁ ∫︁ ∫︁ ∫︁
+ 𝑇𝑚 α𝑁𝑛𝑒 𝑁𝑚
𝑒
𝑑Γ = 𝑁𝑛𝑒 𝑞𝑉 𝑑𝑆 + 𝑁𝑛𝑒 𝑓 𝑑Γ + α𝑁𝑛𝑒 𝑇𝑎 (𝑥)𝑑Γ, (2.5)
Γ3 𝑆 Γ2 Γ3

и уравнения равновесия
∫︁ ∫︁ ∫︁
𝑝1 𝑒
𝑁𝑛,𝑖 𝐶𝑖𝑗𝑘𝑙 ε𝑘𝑙 𝑑𝑆 + 𝑝2 𝑒
𝑁𝑛,𝑖 φ(𝑥, 𝑥′ )𝐶𝑖𝑗𝑘𝑙 ε𝑘𝑙 𝑑𝑆 ′ (𝑥)𝑑𝑆 =
𝑆 𝑆 𝑆 ′ (𝑥)∩𝑆
∫︁ ∫︁ ∫︁
= 𝑝1 𝑒
𝑁𝑛,𝑖 𝐶𝑖𝑗𝑘𝑙 α𝑘𝑙 ∆𝑇 𝑑𝑆 + 𝑝2 𝑒
𝑁𝑛,𝑖 φ(𝑥, 𝑥′ )𝐶𝑖𝑗𝑘𝑙 α𝑘𝑙 ∆𝑇 𝑑𝑆 ′ (𝑥)𝑑𝑆+
𝑆 𝑆 𝑆 ′ (𝑥)∩𝑆
∫︁ ∫︁
+ 𝑁𝑛𝑒 𝑏𝑗 𝑑𝑆 + 𝑁𝑛𝑒 𝑝𝑗 𝑑Γ, (2.6)
𝑆 Γ5

′
где 𝑖,𝑗,𝑘,𝑙 = 1, 2; 𝑚, 𝑛 ∈ 𝐼 𝑒 ; 𝑚′ ∈ 𝐼 𝑒 ; 𝑒 ∈ 𝑆ℎ ; 𝑒′ ∈ 𝑆ℎ′ ; 𝑆ℎ′ — аппроксимированная
зона нелокального влияния, детали аппроксимации которой рассмотрим далее
в следующей главе.
Введём понятия вектора 𝐸 𝑛 — единичный вектор размерности 𝑀 , где
𝑀 — количество узлов в сетке 𝑆ℎ . Переобозначим слагаемые уравнения тепло­
проводности (2.5) при помощи символов
∫︁
𝐿
̂︀ = λ𝑖𝑗 𝑁 𝑒 𝑁 𝑒 𝐸 𝑛 ⊗ 𝐸 𝑚 𝑑𝑆,
K 𝑇 𝑛,𝑖 𝑚,𝑗

∫︁ ∫︁ 𝑆
̂︀ 𝑁 𝐿 =
K 𝑒
𝑁𝑛,𝑖 φ(𝑥, 𝑥′ )λ𝑖𝑗 𝑁𝑚
𝑒 ′
′
′ ,𝑗 𝐸 𝑛 ⊗ 𝐸 𝑚′ 𝑑𝑆 (𝑥)𝑑𝑆,
𝑇
𝑆 𝑆 ′ (𝑥′ )∩𝑆
∫︁
̂︀ α
K = α𝑁𝑛𝑒 𝑁𝑚
𝑒
𝐸 𝑛 ⊗ 𝐸 𝑚 𝑑Γ,
𝑇
Γ3
 30
∫︁
Q= 𝑁𝑛𝑒 𝑞𝑉 𝐸 𝑛 𝑑𝑆,
𝑆
∫︁
F= 𝑁𝑛𝑒 𝑓 𝐸 𝑛 𝑑Γ,
∫︁ Γ2
Tα = α𝑁𝑛𝑒 𝑇𝑎 (𝑥)𝐸 𝑛 𝑑Γ.
Γ3

Аналогичную процедуру сделаем и для уравнения равновесия (2.6)

∫︁
𝐿
̂︀ = 𝑁 𝑒 𝐶𝑖𝑗𝑘𝑙 ε𝑘𝑙 𝐸 𝑛 ⊗ 𝐸 𝑚 𝑑𝑆,
K 𝐸 𝑛,𝑖

∫︁ ∫︁ 𝑆
̂︀ 𝑁 𝐿 =
K 𝑒
𝑁𝑛,𝑖 φ(𝑥, 𝑥′ )𝐶𝑖𝑗𝑘𝑙 ε𝑘𝑙 𝑑𝑆 ′ (𝑥)𝐸 𝑛 ⊗ 𝐸 𝑚′ 𝑑𝑆,
𝐸
𝑆 𝑆 ′ (𝑥)∩𝑆
∫︁
𝐿 𝑒
E
̂︀ = 𝑁𝑛,𝑖 𝐶𝑖𝑗𝑘𝑙 α𝑘𝑙 ∆𝑇 𝐸 𝑛 𝑑𝑆,
∫︁ ∫︁𝑆
̂︀ 𝑁 𝐿 =
E 𝑒
𝑁𝑛,𝑖 φ(𝑥, 𝑥′ )𝐶𝑖𝑗𝑘𝑙 α𝑘𝑙 ∆𝑇 𝐸 𝑛 𝑑𝑆 ′ (𝑥)𝑑𝑆,
𝑆 𝑆 ′ (𝑥)∩𝑆
∫︁
B
̂︀ = 𝑁𝑛𝑒 𝑏𝑗 𝐸 𝑛 𝑑𝑆,
∫︁𝑆
P
̂︀ = 𝑁𝑛𝑒 𝑝𝑗 𝐸 𝑛 𝑑Γ.
Γ5

Тогда после интегрирования, о котором пойдёт речь в следующем разделе, ито­

говые системы можно записать в матрично-векторном виде
𝐿 𝑁𝐿 α
(︁ )︁
𝑝1 K𝑇 + 𝑝2 K𝑇 + K𝑇 · T = Q + F + Tα ,
̂︀ ̂︀ ̂︀ (2.7)
𝐿 𝑁𝐿
̂︀ 𝐿 + 𝑝2 E
̂︀ 𝑁 𝐿 + B
(︁ )︁
𝑝1 K𝐸 + 𝑝2 K𝐸 · U
̂︀ ̂︀ ̂︀ = 𝑝1 E ̂︀ + P.
̂︀ (2.8)

̂︀ 𝐿 и K
Здесь K ̂︀ 𝑁 𝐿 — матрицы локальной и нелокальной теплопроводности;
𝑇 𝑇
̂︀ α — матрица теплообмена; T — вектор искомых узловых значений температу­
K 𝑇

ры; Q и F — векторы дискретизированных внутренних и внешних источников и

̂︀ 𝐿 и K
стоков теплоты; Tα — вектор дискретизированного теплообмена; K ̂︀ 𝑁 𝐿 —
𝐸 𝐸

матрицы локальной и нелокальной жёсткости; U

̂︀ — вектор искомых узловых
 31

перемещений; B
̂︀ и P
̂︀ — векторы дискретизированных плотностей объёмных и
̂︀ 𝐿 и E
поверхностных сил; E ̂︀ 𝑁 𝐿 — векторы локального и нелокального темпера­
̂︀ 𝐿 и K
турного линейного расширения. В силу того, что матрицы K ̂︀ 𝑁 𝐿 имеют
𝐸 𝐸

блочную структуру, с размером блока 2 × 2, для удобства дальнейшего изложе­

ния будем представлять их в виде аналогов (по количеству индексов) тензоров
четвёртого ранга, где первые два индекса обозначают строку и столбец с указа­
нием блока, а вторые — строку и столбец внутри блока. Аналогично представим
векторы U,
̂︀ B,
̂︀ P, ̂︀ 𝐿 и E
̂︀ E ̂︀ 𝑁 𝐿 в виде тензоров второго ранга, где первый индекс

соответствует номеру узла, а второй — номеру координатной компоненты.

2.3. Ассемблирование систем уравнений

Рассмотрим более подробно вопрос ассемблирования систем уравнений

(2.7) и (2.8), которые получены после интегрирования систем (2.5) и (2.6). Для
удобства расмотрим каждое слагаемое отдельно. Но прежде, чем это сделать
̃︀ 𝑒1 𝑒2 и жёсткости K
введём определения блоков матрицы теплопроводности K ̂︀ 𝑒1 𝑒2 ,
𝑛𝑚 𝑛𝑚

стоящих в 𝑛-ой строке и 𝑚-ом столбце соответствующих матриц,

̃︀ 𝑒1 𝑒2 (𝑥, 𝑦) = λ𝑖𝑗 𝑁 𝑒1 (𝑥)𝑁 𝑒1 (𝑦)𝐸 𝑛 ⊗ 𝐸 𝑚 ,

K (2.9)
𝑛𝑚 𝑛,𝑖 𝑚,𝑗

̂︀ 𝑒1 𝑒2 (𝑥, 𝑦) = 𝐶𝑖𝑗𝑘𝑙 𝑁 𝑒1 (𝑥)𝑁 𝑒2 (𝑦)𝐸 𝑛 ⊗ 𝐸 𝑚 ⊗ 𝑒𝑖 ⊗ 𝑒𝑗 ,

K (2.10)
𝑛𝑚 𝑛,𝑘 𝑚,𝑙

где 𝑖,𝑗,𝑘,𝑙 = 1,2; 𝑛,𝑚 = 1,𝑀 . Далее для общности записи будем использовать
блок K𝑒𝑛𝑚
1 𝑒2
, который будет играть роль блока матрицы теплопроводности или
жёсткости в зависимости от контекста.
Рассмотрим матричные слагаемые с множителем 𝑝1 , которые достаточно
легко могут быть аппроксимированы классической конечно-элементной про­
цедурой [72, 138]. После её применения ассемблированную матрицу запишем
следующим образом

̂︀ 𝐿 =
∑︁ ∑︁ ∑︁
K ℱ 𝑤𝑞 K𝑒𝑒 𝑒
𝑛𝑚 (𝑥𝑞 , 𝑥𝑞 )𝐽𝑞 . (2.11)
𝑒∈𝑆ℎ 𝑛,𝑚∈𝐼 𝑒 𝑞∈𝑄𝑒
 32

Здесь 𝑤𝑞 — весовой множитель в квадратурном узле 𝑞; 𝐽𝑞𝑒 — аппроксими­

рованный якобиан в квадратурном узле 𝑞 на элементе 𝑒; 𝑥𝑞 — коордианата
квадратурного узла под номером 𝑞; 𝑄𝑒 — набор номеров квадратурных узлов
на элементе 𝑒.
Аппроксимацию интегральных слагаемых, стоящих у множителей 𝑝2 урав­
нений (2.5) и (2.6), следует начать с аппроксимации зоны нелокального влияния
𝑆ℎ′ , для которой необходимо вначале аппроксимировать внешние интегралы, где
приходим к промежуточным выражениям следующего вида
∫︁
𝑁𝐿 ∑︁ ∑︁ ∑︁
𝑒′
K
̂︀
𝑇 =
𝑒 𝑒
𝑤𝑞 𝑁𝑛,𝑖 (𝑥𝑞 )𝐽𝑞 φ(𝑥, 𝑥′ )λ𝑖𝑗 𝑁𝑚 ′
′ ,𝑗 𝑑𝑆 (𝑥)𝐸 𝑛 ⊗ 𝐸 𝑚′ ,

𝑒∈𝑆ℎ 𝑛∈𝐼 𝑒 𝑞∈𝑄𝑒

𝑆 ′ (𝑥′ )∩𝑆
∫︁
̂︀ 𝑁 𝐿
∑︁ ∑︁ ∑︁
K 𝐸 = 𝑒
𝑤𝑞 𝑁𝑛,𝑖 (𝑥𝑞 )𝐽𝑞𝑒 φ(𝑥, 𝑥′ )𝐶𝑖𝑗𝑘𝑙 ε𝑘𝑙 𝑑𝑆 ′ (𝑥)𝐸 𝑛 ⊗ 𝐸 𝑚′ .
𝑒∈𝑆ℎ 𝑛∈𝐼 𝑒 𝑞∈𝑄𝑒
𝑆 ′ (𝑥′ )∩𝑆

Далее в каждом квадратурном узле 𝑥𝑞 необходимо аппроксимировать область

нелокального влияния 𝑆ℎ𝑞 [120], которую можно представить в виде множества
элементов, квадратурные узлы которых хотя бы частично попали в область
𝑆 ′ (𝑥𝑞 ). На Рис. 2.1 крестом указан узел относительно которого проводит­
ся аппроксимация, область нелокального влияния ограничена окружностью,
точками отмечены квадратурные узлы элементов, а серым цветом выделены
элементы, которые были учтены в аппроксимированной области нелокально­
го влияния 𝑆ℎ𝑞 . Такой способ аппроксимации будем называть квадратурной
аппроксимацией. Тогда ассемблирование матрицы соответствующей нелокаль­
ному слагаемому запишем следующим образом

̂︀ 𝑁 𝐿 = ′ ′
∑︁ ∑︁ ∑︁ ∑︁ ∑︁ ∑︁
K ℱ 𝑤𝑞 𝐽𝑞𝑒 𝑤𝑞′ φ(𝑥𝑞 , 𝑥𝑞′ )K𝑒𝑒 𝑒
𝑛𝑚′ (𝑥𝑞 , 𝑥𝑞 ′ )𝐽𝑞 ′ . (2.12)
𝑒∈𝑆ℎ 𝑛∈𝐼 𝑒 𝑞∈𝑄𝑒 𝑒′ ∈𝑆ℎ𝑞 𝑚′ ∈𝐼 𝑒′ 𝑞 ′ ∈𝑄𝑒′

Ассемблирование остальных слагаемых уравнения теплопроводности (2.7)

происходит без каких-либо особенностей, поэтому просто выпишем их без по­
 33

Рис. 2.1. Квадратурная аппроксимация области нелокального влияния

дробного разъяснения деталей

α ∑︁ ∑︁ ∑︁
K𝑇 =
̂︀ 𝑤𝑞 α𝑁𝑛𝑒 (𝑥𝑞 )𝑁𝑚
𝑒
(𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 ⊗ 𝐸 𝑚 , (2.13)
𝑒∈Γℎ 𝑛,𝑚∈𝐼 𝑒 𝑞∈𝑄𝑒
∑︁ ∑︁ ∑︁
Tα = 𝑤𝑞 α𝑁𝑛𝑒 (𝑥𝑞 )𝑇α (𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 , (2.14)
𝑒∈Γℎ 𝑛∈𝐼 𝑒 𝑞∈𝑄𝑒
∑︁ ∑︁ ∑︁
Q= 𝑤𝑞 𝑞𝑉 (𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 , (2.15)
𝑒∈𝑆ℎ 𝑛∈𝐼 𝑒 𝑞∈𝑄𝑒
∑︁ ∑︁ ∑︁
F= 𝑤𝑞 𝑓 (𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 . (2.16)
𝑒∈Γℎ 𝑛∈𝐼 𝑒 𝑞∈𝑄𝑒

Аналогично запишем ассемблирование остальных слагаемых и для уравнения

равновесия (2.8)
∑︁ ∑︁ ∑︁
B
̂︀ = 𝑤𝑞 𝑏(𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 , (2.17)
𝑒∈𝑆ℎ 𝑛∈𝐼 𝑒 𝑞∈𝑄𝑒
∑︁ ∑︁ ∑︁
P
̂︀ = 𝑤𝑞 𝑝(𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 , (2.18)
𝑒∈Γℎ 𝑛∈𝐼 𝑒 𝑞∈𝑄𝑒

̂︀ 𝐿 =
∑︁ ∑︁ ∑︁
E 𝑤𝑞 ∇𝑁𝑛𝑒 (𝑥𝑞 )C
̂︀ · ·α∆𝑇
̂︀ (𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 , (2.19)
𝑒∈𝑆ℎ 𝑛∈𝐼 𝑒 𝑞∈𝑄𝑒

̂︀ 𝑁 𝐿 =
∑︁ ∑︁ ∑︁
E 𝑤𝑞 ∇𝑁𝑛𝑒 (𝑥𝑞 )𝐽𝑞𝑒 ×
𝑒∈𝑆ℎ 𝑛∈𝐼 𝑒 𝑞∈𝑄𝑒
′
∑︁ ∑︁
× ̂︀ · ·α∆𝑇
𝑤𝑞′ φ(𝑥𝑞 , 𝑥𝑞′ )C ̂︀ (𝑥𝑞′ )𝐽𝑞𝑒′ 𝐸 𝑛 . (2.20)
𝑒′ ∈𝑆ℎ𝑞 𝑞 ′ ∈𝑄𝑒′
 34

̂︀ 𝑁 𝐿 был использован метод квадратурной

Отметим, что при ассемблировании E
аппроксимации области нелокального влияния 𝑆 ′ (𝑥), который ранее применял­
ся к матрицам нелокальной телопроводности и жёсткости (2.12).

2.4. Вычисление производных величин

После решения СЛАУ (2.7) и (2.8) на основе полученных сеточных функ­

ций температуры T и перемещения U
̂︀ можем найти их производные величины,

такие как вектор плотности теплового потока 𝑞 и тензор напряжений σ

̂︀ соот­
ветственно. Для этого вычислим градиенты сеточных функций в квадратурных
узлах, пользуясь при этом формулами (2.3) и (2.4). Далее аппроксимируем ин­
тегралы (1.3) и (1.6), для этого снова воспользуемся процедурой квадратурной
аппроксимации области нелокального влияния, после чего получаем формулы
для вычисления вектора плотности теплового потока
⎛ ⎞
𝑒′ ⎠
∑︁ ∑︁
𝑒 𝑒
𝑞 𝑞 = ⎝−𝑝1 λ𝑇𝑚 𝑁𝑚,𝑘 (𝑥𝑞 ) − 𝑝2 𝑤𝑞′ λ𝑇𝑚′ 𝑁𝑚 ′ ,𝑘 (𝑥𝑞 ′ )𝐽𝑞 ′ 𝑒𝑘 (2.21)
𝑒′ ∈𝑆ℎ𝑞 𝑞 ′ ∈𝑄𝑒′

и тензора напряжений
(︃
̂︀ 𝑞 = 𝑝1 𝐶𝑖𝑗𝑘𝑙 (ε𝑘𝑙 (𝑥𝑞 ) − α𝑘𝑙 ∆𝑇𝑞 ) +
σ
)︃
′
∑︁ ∑︁
+ 𝑝2 𝑤𝑞′ 𝐶𝑖𝑗𝑘𝑙 (ε𝑘𝑙 (𝑥𝑞′ ) − α𝑘𝑙 ∆𝑇𝑞′ ) 𝐽𝑞𝑒′ 𝑒𝑘 ⊗ 𝑒𝑙 , (2.22)
𝑒′ ∈𝑆ℎ𝑞 𝑞 ′ ∈𝑄𝑒′

где 𝑞 𝑞 , σ
̂︀ 𝑞 и ∆𝑇𝑞 — значения вектора плотности теплового потока, тензора на­
пряжений и разницы температур в квадратурном узле 𝑞 соответственно.
Для дальнейшего анализа переинтерполируем решения из квадратурных
узлов в регулярные узлы сетки. Для этого, при рассмотрении конкретного уз­
ла сетки, определим ближайшие квадратурные узлы каждого из элементов, в
состав которых входит этот узел, после чего вычислим среднюю величину с
учётом площадей элементов. Другими словами, вначале для каждого элемента
 35

𝑒 необходимо решить задачу минимизации с поиском нужного индекса квадра­

турного узла на элементе 𝑞 𝑒

min𝑒 ρ(𝑥𝑛 , 𝑥𝑞 ) → 𝑞 𝑒 , 𝑛 ∈ 𝑆ℎ .
𝑞∈𝑄

После чего проводим процедуру осреднения с весами, где в качестве весовых

множителей используем площади элементов
∑︁ ∑︁ 𝑞 𝑞𝑒 |𝑆 𝑒 | ∑︁ σ̂︀ 𝑞𝑒 |𝑆 𝑒 |
𝑛 𝑒
|𝑆 | = |𝑆 |, 𝑞𝑛 = 𝑛|
, σ
̂︀ 𝑛 = 𝑛|
,
|𝑆 |𝑆
𝑒∈𝐸 𝑛 𝑛
𝑒∈𝐸 𝑛
𝑒∈𝐸

где 𝑞 𝑛 и σ
̂︀ 𝑛 — значения плотности теплового потока и напряжений в регулярных
узлах сетки соответственно.

2.5. Основные результаты и выводы по главе 2

1. На основе метода конечных элементов, разработана численная схема ап­

проксимации уравнений теплопроводности и равновесия в нелокальных
постановках.
2. Предложен способ квадратурной аппроксимации области нелокального
влияния, суть которого заключена в аппроксимации области нелокаль­
ного влияния относителько каждого квадратурного узла.
3. Представлен алгоритм аппроксимации производных величин, таких как
вектор плотности теплового потока и тензор напряжений, с учётом про­
странственной нелокальности.
 36

Глава 3. Реализация программного комплекса

3.1. Общая структура программного комплекса

В рамках диссертационной работы был реализован конечно-элементный

программный комплекс NonLocFEM [56]. Основная задача комплекса — эф­
фективное решение термомеханических задач в нелокальных постановках на
современных вычислительных системах с использованием технологий парал­
лельных вычислений OpenMP [13] и MPI [12]. Все описанные далее методы
и алгоритмы реализованы в рамках данного комплекса, а именно: аппроксима­
ция области нелокального влияния; параллельные алгоритмы ассемблирования
конечно-элементных матриц; алгоритмы балансировки данных между процесса­
ми и потоками исполнения; интегрирование с использованием нестандартных
базисов конечных элементов; решатели СЛАУ с использованием специально
разработанных предобуславливателей; а также многие другие алгоритмы и ме­
тоды, на которых не будем заострять слишком много внимания.
Глобальная структура программного комплекса включает в себя мате­
матическое ядро и обработчик конфигурационных файлов. Математическое
ядро, в свою очередь, также состоит из нескольких взаимосвязанных библио­
тек, где в качестве основных можно выделить следующие: metamath, parallel,
mesh и solvers. В них находятся необходимые примитивы и алгоритмы для ко­
нечно-элементных расчётов. Обработчик конфигурационных файлов работает
со структурами, представленными в формате JSON [10] и на их основе формиру­
ет запросы для математического ядра, которое проводит необходимые расчёты
и возвращает результаты в форматах, которые можно прочитать популярны­
ми программами для визуализации данных, например Paraview [14]. Помимо
собственных разработок в зависимости комплекса входят две сторонние биб­
лиотеки: библиотека линейной алгебры Eigen [7] и библиотека для работы со
 37

структурами в формате JSON N. Lohmann [11]. Схема взаимосвязи модулей

программы представлена на Рис. 3.1, где зависимый модуль указывает стрел­
кой на модуль от которого он зависит.

Математическое Обработчик
ядро конфигурационных
файлов
Metamath Parallel

Mesh nlohmann/json

Eigen Solvers Configs

NonLocFEM
Рис. 3.1. Структура программы NonLocFEM

Программный комплекс NonLocFEM реализован на языке программиро­

вания C++ [3]. Выбор в пользу этого языка был обоснован его популярностью,
богатой стандартной библиотекой, а главное производительностью итоговых
программ. Помимо этого, язык предоставляет широкий спектр возможностей в
реализации своих идей, особенно если говорить про актуальный на сегодняшний
день стандарт языка C++23, возможности которого повсеместно использова­
ны в программном комплексе. Язык C++ мультипарадигменный, поэтому в
нём существует возможность совмещать объектно-ориентированные и функци­
ональные подходы к программированию. Объектно-ориентированный подход
выражен в виде возможности создания достаточно сложных иерархий классов,
в которых могут быть использованы виртуальные методы. Это в свою очередь
подразумевает позднее связывание кода программы, то есть объекты с разной
логикой обработки тех или иных данных, имеющие при этом единый интерфейс,
 38

могут быть созданы динамически во время выполнения программы. Функци­

ональная парадигма в контексте языка программирования C++ выражена в
виде метапрограммирования шаблонов [19, 26], то есть статична, где часть
вычислений можно вынести на этап компиляции программы, на основе результа­
тов которых генерируется конечный исполняющий файл. Комбинирование двух
парадигм открывает возможность совместить такие, порой несовместимые, ас­
пекты программы, как гибкость исходного кода с его производительностью.
Учитывая специфику вычислительных программ, в разработке программного
комплекса NonLocFEM в большей степени было отдано предпочтение функци­
ональным подходам к программированию.
Наибольшее количество приёмов метапрограммирования было задейство­
вано в библиотеке metamath, за счёт чего она и получила такое название. В этой
библиотеке реализованы различные математические примитивы и функции,
а также представлена адаптированная версия библиотеки символьного диф­
ференцирования на этапе компиляции symdiff, основную концепцию которой
можно найти в монографии Краснова М.М. [36]. Библиотека symdiff содержит
в себе базовые примитивы: константа, переменная; математические операции,
такие как, сложение, вычитание, умножение, деление; ряд математических
функций, включающих в себя экспоненту, логарифм, тригонометрические
функции и многие другие. Благодаря этим примитивам можно строить выра­
жения любой сложности, а также комбинировать эти выражения между собой.
Каждое такое выражение образует уникальный тип данных. Все эти типы дан­
ных объединяет общий интерфейс, предоставляющий возможность вычислить
значение этого выражения в точке и посчитать его производную по заданной
переменной, при этом вычисление производной порождает новый тип данных,
который регистрируется на этапе компиляции программы и соотвествует вы­
ражению, являющимся производной исходного выражения. При этом были
 39

предприняты меры по оптимизации конечных выражений, сокращающие ко­

личество операций, которые необходимы при вычислении значения в точке.
На основе библиотеки symdiff, в рамках библиотеки metamath, была
построена библиотека конечных элементов finite_elements, в которой symdiff
использована при описании базисных функций форм элементов и вычислении
их производных компилятором. Такой подход значительно упрощает процесс
добавления новых элементов, а также сокращает время отладки программы,
так как большая часть ошибок, как правило, происходит именно при ручном
дифференцировании базисов элементов. Помимо этого такой подход ещё со­
кращает исходный код программы, так как при добавлении нового конечного
элемента прикладному программисту необходимо всего лишь описать базис но­
вого элемента в терминах библиотеки symdiff и записать координаты его узлов
в локальной системе координат элемента. Для некоторых семейств элементов
такой процесс можно автоматизировать, например, для семейства лагранже­
вых элементов, где базисы элементов построены по определённому алгоритму.
Таким образом, можно указать порядок элемента, после чего компилятор сге­
нерирует необходимый базис в виде набора выражений, которые в дальнейшем
могут быть им же и продифференцированы. Помимо дифференцирования бази­
сов, в этой библиотеке также представлен функционал, связанный с процедурой
интегрирования. Интегрирование, в отличие от дифференцирования, здесь ре­
ализовано численно. Библиотека содержит в себе наборы квадратур разного
порядка, при помощи которых выполняется процедура интегрирования. Все эле­
менты и квадратуры, описанные в рамках данной библиотеки, имеют единый
интерфейс, поэтому дальнейшие конечно-элементные алгоритмы могут быть
записаны в обобщённой форме.
Библиотека mesh предназначена для работы с конечно-элементными
сетками. Основным классом данной библиотеки является класс хранилище,
объекты которого могут читать файлы с сетками и представлять их в виде,
 40

с которым взаимодествуют алгоритмы программы, в том числе и алгоритмы

модуля solvers. Схема хранения подразумевает, что элементы образованы пу­
тём перечисления номеров узлов, которые им принадлежат, а также ссылкой
на объект библиотеки finite_elements, в котором определены функции формы
и квадратурные узлы в локальной системе координат элемента. Также в схеме
хранения участвуют координаты узлов сетки и именованные группы элемен­
тов, образующих подобласти для определения границ, на которых далее можно
задать разные граничные условия. Также было принято решение о том, что
в случае использования распределённых вычислений, при помощи библиотеки
MPI, сетка представлена целиком на каждом отдельном процессе выполнения
программы. Такое решение связано с желанием не усложнять схему хранения,
ведь для задач в нелокальных постановках, как правило, используются доста­
точно грубые сетки, содержащие малое количество элементов. Однако, это не
касается алгоритмов аппроксимации области нелокального влияния, в котором
задействован алгоритм поиска ближайших соседей. В результате алгоритма по­
иска ближайших соседей хранятся только те данные, которые необходимы для
ассемблирования куска матрицы, обрабатываемой конкретным процессом.
Помимо класса хранилища, в библиотеке mesh также реализованы ал­
горитмы для работы с конечно-элеметными сетками. Здесь реализованы ал­
горитмы вычисления квадратурных узлов в глобальной системе координат,
вычисления в них якобианов и производных функций форм относительно гло­
бальных переменных. Здесь реализован шаблонный параллельный алгоритм
обхода по сетке, на базе которого построены алгоритмы балансировки данных,
алгоритм перенумерации узлов Катхилла — Макки [78], алгоритм аппроксима­
ции области нелокального влияния и алгоритмы формирования и заполнения
портрета конечно-элеметных матриц из модуля solvers. Также в библиотеке
mesh реализован функционал сохранения результатов расчётов.
 41

В библиотеке solvers реализованы алгоритмы ассемблирования матриц и

правых частей, а также решатели СЛАУ. Алгоритмы ассемблирования вклю­
чают в себя алгоритмы формирования портрета матриц и его заполнения, при
этом алгоритмы работают построчно, что обеспечивает хорошую масштабиру­
емость на многопроцессорных системах. Здесь же представлены обобщённые
решатели задач теплопроводности и равновесия, на вход которым подаются
расчётная сетка, параметры материала и граничные условия. После чего дан­
ные решатели возвращают решения, содержащие искомые величины, такие как
температуру и перемещения, и производные от них тепловые потоки, дефор­
мации и напряжения.
Продолжая обзор библиотек программы следует поговорить о библиотеке
parallel. Во многом эта библиотека является «обёрткой» над библиотеками па­
раллельного программирования OpenMP и MPI. Она служит двум основным
принципам: во-первых — осовременить и обобщить интерфейсы давно устов­
шихся функций из ранее упомянутых библиотек, которые в современном C++
выглядят весьма громоздко, а во-вторых — собрать все параллельные вызовы
в одном месте, что даёт возможность легко отключать параллелилизм без до­
полнительных сложностей связанных с компиляцией программы. Кроме того
эта библиотека содержит в себе основы для алгоритмов балансировки данных,
которые затем используются в более общих алгоритмах библиотеки mesh, и
удобные примитивы для работы с параллельным кодом.
Для связи структур описанных в формате JSON с математическим яд­
ром была разработана библиотека configs. Она содержит в себе примитивы,
которые вычисляются на основе JSON и интерпретируются конечной програм­
мой в терминах представленных в математическом ядре программы. В случае
несоответствия структуры описанной в JSON с ожидаемой библиотека генери­
рует сообщения об ошибках с указанием места в конфигурационном файле,
 42

где возникла эта ошибка. Пример конфигурационного файла с описанием его

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

3.2. Параллельный алгоритм ассемблирования матриц

Задачи в нелокальных постановках обладают достаточно большой вы­

числительной сложностью. Матрицы, получаемые после конечно-элементной
аппроксимации, значительно более плотные по сравнению с их классически­
ми аналогами, а также требуют огромных вычислительных ресурсов для их
ассемблирования [104]. Поэтому возникает спрос на использование всех воз­
можностей современных компьютеров, а именно параллельные вычисления на
машинах с общей и распределённой памятью. Но для того, чтобы полностью
задействовать все вычислительные ресурсы, необходимо разработать алгоритм
пригодный для распараллеливания.
Для начала упростим аппроксимацию области нелокального влияния
𝑆 ′ (𝑥), так как рассматриваемый ранее способ квадратурной аппроксимации,
описанный формулой (2.12) и представленный на Рис. 2.1, крайне неудобен для
практической реализации. Поэтому возникает идея упростить его и проводить
аппроксимацию не относительно квадратурных узлов сетки, а относительно цен­
тров элементов. На Рис. 3.2 крестом обозначен центр элемента относительно
которого проводится аппроксимация, область нелокального влияния ограни­
чена окружностью, точками обозначены центры элементов, а серым цветом
выделены элементы образующие аппроксимированную область нелокального
 43

влияния 𝑆ℎ𝑒 . Такой способ аппроксимации будем называть элеметной аппрок­

симацией. Тогда в алгоритме ассемблирования нелокальной матрицы (2.12)
сможем поменять местами знаки суммирования, после чего сможем преобра­
зовать его к следующему виду

̂︀ 𝑁 𝐿 = ′ ′
∑︁ ∑︁ ∑︁ ∑︁ ∑︁ ∑︁
K ℱ 𝑤𝑞 𝐽𝑞𝑒 𝑤𝑞′ φ(𝑥𝑞 , 𝑥𝑞′ )K𝑒𝑒 𝑒
𝑛𝑚 (𝑥𝑞 , 𝑥𝑞 ′ )𝐽𝑞 ′ . (3.1)
𝑒∈𝑆ℎ 𝑛∈𝐼 𝑒 𝑒′ ∈𝑆ℎ𝑒 𝑚∈𝐼 𝑒′ 𝑞∈𝑄𝑒 𝑞 ′ ∈𝑄𝑒′

Также можем упростить и алгоритм вычисления нелокального температурного

линейного расширения (2.20)

̂︀ 𝑁 𝐿 = ′
∑︁ ∑︁ ∑︁ ∑︁ ∑︁
E 𝑤𝑞 ∇𝑁𝑛𝑒 (𝑥𝑞 )𝐽𝑞𝑒 ̂︀ · ·α∆𝑇
𝑤𝑞′ φ(𝑥𝑞 , 𝑥𝑞′ )C ̂︀ (𝑥𝑞′ )𝐽𝑞𝑒′ 𝐸 𝑛 .
𝑒∈𝑆ℎ 𝑛∈𝐼 𝑒 𝑒′ ∈𝑆ℎ𝑒 𝑞∈𝑄𝑒 𝑞 ′ ∈𝑄𝑒′

(3.2)

Рис. 3.2. Элементная аппроксимация области нелокального влияния

Такой подход позволяет отделить алгоритм обхода по сетке, ответственно­

го за формирование портрета матрицы, от алгоритма интегрирования, однако,
он обладает дефектом, который заключается в том, что не все квадратурные
узлы, попадающие в область 𝑆 ′ (𝑥𝑞 ), центр которой находится в некотором квад­
ратурном узле под номером 𝑞, попадают под покрытие аппроксимированной
области нелокального влияния 𝑆ℎ𝑒 . Это может приводить к нарушениям балан­
са, что в свою очередь приводит к менее точному решению и даже осциляциям.
 44

Для решения данной проблемы, радиус поиска соседних элементов нужно брать
больше радиуса нелокальности 𝑟, например, на величину максимального рассто­
яния между центрами двух смежных элементов, где под смежными элементами
подразумеваем элементы обладающие хотя бы одним общим узлом. Таким об­
разом, все необходимые квадратурные узлы будут учтены в расчёте.
После разделения алгоритма обхода сетки и алгоритма интегрирования,
можем изменить первый таким образом, чтобы сделать его пригодным для
параллельных вычислений. Главной проблемой алгоритма (3.1) остаётся за­
висимость номера узла сетки от номера текущего элемента из-за чего при
использовании параллельных вычислений существует вероятность возникнове­
ния гонки данных, что, в зависимости от подхода к распараллеливанию, может
приводить к неправильному решению задачи, или частым барьерным синхро­
низациям, которые в свою очередь снижают эффективность использования
параллельных вычислений. Поэтому возникает идея изменить порядок сумми­
рования таким образом, чтобы такой зависимости не было. Для этого определим
для каждого узла сетки 𝑛 ∈ 𝑆ℎ множество элементов 𝐸 𝑛 , которым он принад­
лежит и изменим порядок суммирования так, чтобы под первым знаком суммы
были номера узлов, а под вторым номера элементов которым он принадлежит

̂︀ 𝑁 𝐿 = ′ ′
∑︁ ∑︁ ∑︁ ∑︁ ∑︁ ∑︁
K ℱ 𝑤𝑞 𝐽𝑞𝑒 𝑤𝑞′ φ(𝑥𝑞 , 𝑥𝑞′ )K𝑒𝑒 𝑒
𝑛𝑚 (𝑥𝑞 , 𝑥𝑞 ′ )𝐽𝑞 ′ . (3.3)
𝑛∈𝑆ℎ 𝑒∈𝐸 𝑛 𝑒′ ∈𝑆ℎ𝑒 𝑚∈𝐼 𝑒′ 𝑞∈𝑄𝑒 𝑞 ′ ∈𝑄𝑒′

Такой алгоритм сборки матрицы является построчным, соответственно каждую

строку матрицы можно собирать независимо в своём исполняемом потоке, а
также распределить вычисление строк между вычислительными узлами.
Полученный алгоритм (3.3) пригоден для параллельных вычислений на
машинах с общей и распределённой памятью, однако, возникает проблема ба­
лансировки данных и объёма вычислений. При решении этих проблем стоит
начинать с проблемы балансировки данных, так как из-за неё есть вероятность
возникновения ситуации, когда задача не может быть решена, в силу того, что
 45

на одном из вычислительных узлов может не хватить оперативной памяти, в то

время как при балансировке данных такой проблемы можно было бы избежать.
При балансировке данных будем исходить из гипотезы, что на каждом
вычислительном узле 𝑝 ∈ 𝑃 установлено одинаковое количество оперативной
памяти и данные нужно распределить равномерно между всеми вычислитель­
ными узлами 𝑃 . Для этого необходимо найти общее число элементов матрицы
𝑀 , затем найти среднее 𝑀𝑚 = 𝑀/|𝑃 | и распределить строки матрицы между
вычислительными узлами таким образом, чтобы в каждой группе строк количе­
ство элементов матрицы было приблизительно равным 𝑀𝑝 ≈ 𝑀𝑚 . На практике
удобнее всего брать группы последовательных строк. Причём при балансиров­
ке нет необходимости формировать полный портрет матрицы, это также можно
делать построчно, что гораздо эффективнее и легко реализовать. Также при ба­
лансировке необходимо учесть симметрию полученной матрицы, так как из-за
достаточно больших объёмов выгодно хранить лишь её половину.
После балансировки данных между вычислительными узлами, можем
также провести балансировку объёмов вычислений между вычислительными
потоками. Для этого необходимо проделать ту же процедуру, что и при ба­
лансировке данных, но осреднять не по количеству элементов матрицы, а по
количеству вызовов функции интегрирования. В некоторых ситуациях такая
балансировка может быть полезна и между вычислительными узлами, напри­
мер, когда сборка матрицы занимает гораздо больше времени, чем решение
итоговой СЛАУ.
Вместе с параллельным алгоритмом ассемблирования нелокальных мат­
риц (3.3), можем также выписать параллельный алгоритм ассемблирования
локальных матриц (2.11)

̂︀ 𝐿 =
∑︁ ∑︁ ∑︁ ∑︁
K ℱ 𝑤𝑞 K𝑒𝑒 𝑒
𝑛𝑚 (𝑥𝑞 , 𝑥𝑞 )𝐽𝑞 ,
𝑛∈𝑆ℎ 𝑒∈𝐸 𝑛 𝑚∈𝐼 𝑒 𝑞∈𝑄𝑒
 46

Аналогично можно поступить и с ассемблированием матрицы теплообмена

(2.13), а также векторами в правой части (2.14) — (2.16) уравнения теплопро­
водности (2.7)
∑︁ ∑︁ ∑︁
Q= 𝑤𝑞 𝑞𝑉 (𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 ,
𝑛∈𝑆ℎ 𝑒∈𝐸 𝑛 𝑞∈𝑄𝑒

∑︁ ∑︁ ∑︁
F= 𝑤𝑞 𝑓 (𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 ,
𝑛∈Γℎ 𝑒∈𝐸 𝑛 𝑞∈𝑄𝑒

̂︀ α =
∑︁ ∑︁ ∑︁ ∑︁
K 𝑇 𝑤𝑞 α𝑁𝑛𝑒 (𝑥𝑞 )𝑁𝑚
𝑒
(𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 ⊗ 𝐸 𝑚 ,
𝑛∈Γℎ 𝑒∈𝐸 𝑛 𝑚∈𝐼 𝑒 𝑞∈𝑄𝑒

∑︁ ∑︁ ∑︁
Tα = 𝑤𝑞 α𝑁𝑛𝑒 (𝑥𝑞 )𝑇α (𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 .
𝑛∈Γℎ 𝑒∈𝐸 𝑛 𝑞∈𝑄𝑒

Проделаем те же выкладки и для векторов правой части (2.17) — (2.19) и (3.2)

уравнения равновесия (2.8)
∑︁ ∑︁ ∑︁
B
̂︀ = 𝑤𝑞 𝑏(𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 ,
𝑛∈𝑆ℎ 𝑒∈𝐸 𝑛 𝑞∈𝑄𝑒

∑︁ ∑︁ ∑︁
P
̂︀ = 𝑤𝑞 𝑝(𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 ,
𝑒∈Γℎ 𝑛∈𝐼 𝑒 𝑞∈𝑄𝑒

̂︀ 𝐿 =
∑︁ ∑︁ ∑︁
E 𝑤𝑞 ∇𝑁𝑛𝑒 (𝑥𝑞 ) · C
̂︀ · ·α∆𝑇
̂︀ (𝑥𝑞 )𝐽𝑞𝑒 𝐸 𝑛 ,
𝑛∈𝑆ℎ 𝑒∈𝐸 𝑛 𝑞∈𝑄𝑒

̂︀ 𝑁 𝐿 = ′
∑︁ ∑︁ ∑︁ ∑︁ ∑︁
E 𝑤𝑞 ∇𝑁𝑛𝑒 (𝑥𝑞 )𝐽𝑞𝑒 · ̂︀ · ·α∆𝑇
𝑤𝑞′ φ(𝑥𝑞 , 𝑥𝑞′ )C ̂︀ (𝑥𝑞′ )𝐽𝑞𝑒′ 𝐸 𝑛 ,
𝑛∈𝑆ℎ 𝑒∈𝐸 𝑛 𝑒′ ∈𝑆ℎ𝑒 𝑞∈𝑄𝑒 𝑞 ′ ∈𝑄𝑒′

Отметим, что при ассемблировании векторов на каждом вычислительном уз­

ле можем также ограничиться только теми строками, которые были получены
при балансировке.
 47

3.3. Алгоритм аппроксимации области нелокального влияния

При аппроксимации области нелокального влияния необходимо использо­

вать метод поиска ближайших соседей, однако, важно выбрать оптимальный
для рассматриваемых задач алгоритм. Наивный алгоритм линейного поиска
имеет квадратичную сложность 𝑂(𝑁 2 ), в связи с чем на достаточно подроб­
ных сетках время поиска может быть весьма существенным. Поэтому построим
алгоритм поиска в основе которого лежит k-d дерево [73].
Суть метода на основе k-d дерева заключается в разделении простран­
ства занимаемого телом на равномерные ячейки, после разбиения на которые
необходимо составить списки узлов соответствующие ячейкам, в которых они
оказались. Размер ячеек следует брать равным радиусу поиска, тогда при поис­
ке ближайших соседей поиск можно ограничить до ячейки в которой находится
этот узел и смежных ей ячейкам. В качестве алгоритма поиска внутри яче­
ек можно использовать обычный линейный поиск, тогда сложность такого
алгоритма можно оценить как 𝑂(𝑁 log 𝑁 ), что на подробных сетках замет­
но быстрее линейного поиска. Пример разбиения представлен на Рис. 3.3, где
крестом указан рассматриваемый узел, тёмно-серым цветом выделена ячейка,
которой принадлежит этот узел, а светло-серым смежные ей ячейки, область
нелокального влияния ограничена окружностью.

3.4. Оптимизация базисных функций конечных элементов

Помимо эффективных алгоритмов сборки матрицы жёсткости, также воз­

никает потребность в эффективном решении итоговых СЛАУ. Прямые методы,
применённые к разреженным матрицам, как правило требуют значительных
затрат оперативной памяти, поэтому возникает спрос на использование ите­
рационных методов решения [31]. Однако итерационные методы могут иметь
слишком медленную сходимость, которая в первую очередь обусловлена самой
 48

Рис. 3.3. K-d дерево для поиска ближайших соседей

системой. Таким образом возникает идея уменьшить число обусловленности си­

стемы и как правило для этого используют разного рода предобуславливатели
и адаптируют сетку таким образом, чтобы элементы были ближе к своим ис­
ходным формам [108, 119]. Также для ускорения сходимости можно подобрать
начальные данные, чтобы они были как можно точнее к искомому решению.
Но в случае, когда для вычислений используют элементы высшего порядка,
появляется дополнительная возможность оптимизировать базис элемента под
конкретную задачу.
Обычно использование элементов высшего порядка может быть связано с
желанием более точно аппроксимировать искомые величины или производные
от этих величин, такие как тепловые потоки или деформации [30, 80]. Также они
могут понадобиться для решения специфических задач, где могут встречаться
производные высшего порядка. В задачах не обладающих такой спецификой
большой популярностью обладают квадратичные серендиповые элементы, так
как они позволяют достаточно точно аппроксимировать градиенты искомых
функций и при этом, как показано на Рис. 3.4, не имеют внутренних узлов,
что заметно упрощает расчёты.
 49

7 6 5

ξ
8 4
O

1 2 3

Рис. 3.4. Квадратичный серендиповый элемент в локальной системе координат

Oξη

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

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

1
𝑁𝑖 = (1 + ξ𝑖 ξ)(1 + η𝑖 η)((9𝑠 − 1)(1 − ξ𝑖 ξ − η𝑖 η) + (9𝑠 + 3)ξ𝑖 ξη𝑖 η),
16
𝑖 = 1, 3, 5, 7; ξ𝑖 , η𝑖 = ±1,
1
𝑁𝑖 = (1 − ξ2 )(1 + η𝑖 η)((5 − 9𝑠) + (9𝑠 + 3)η𝑖 η), 𝑖 = 2, 6; η𝑖 = ±1,
16
1
𝑁𝑖 = (1 − η2 )(1 + ξ𝑖 ξ)((5 − 9𝑠) + (9𝑠 + 3)ξ𝑖 ξ), 𝑖 = 4, 8; ξ𝑖 = ±1.
16

Выбор параметра 𝑠 был основан на следующих предположениях

∫︁1 ∫︁1
𝑁𝑖 (ξ, η)𝑑ξ𝑑η = 𝑠, 𝑖 = 1, 3, 5, 7,
−1 −1
∫︁ ∫︁1
1

𝑁𝑖 (ξ, η)𝑑ξ𝑑η = 1 − 𝑠, 𝑖 = 2, 4, 6, 8.
−1 −1
 50

Таким образом, классический базис можно получить при 𝑠 = −1/3.

Теперь можем перейти к задаче минимизации числа обусловленности мат­
риц теплопроводности K
̂︀ 𝑇 и жёсткости K
̂︀ 𝐸 , которые для удобства дальнейшего

изложения будем обозначать одной буквой K.

̂︀ Для этого введём понятие числа

обусловленности, как квадратный корень отношения максимального по моду­

лю собственного числа λmax к минимальному по модулю собственному числу
λmin матрицы K
̂︀
√︃
|λmax |
cond K
̂︀ = . (3.4)
|λmin |
Воспользуемся гипотезой, что минимальное по модулю собственное число
λ𝑚𝑖𝑛 слабо зависит от параметров модели и дополнительного параметра базиса.
Так как след матрицы равен сумме её собственных значений, а сама матрица
симметричная и положительно определённая, то задача минимизации числа
обусловленности эквивалентна задаче минимизации следа матрицы
𝑘 ∑︁
∑︁ ∫︁ ∑︁ 2
𝑒
min tr K
̂︀ = min 𝑐𝑗 𝑁𝑖,𝑗 𝑑𝑆 𝑒 ,
𝑠 𝑠
𝑒∈𝑆ℎ 𝑆 𝑒 𝑖=0 𝑗=0

где 𝑘 — количество функций форм элемента, 𝑐𝑗 — постоянные коэффициенты,

зависящие от свойств материала. Но, как легко заметить, результат задачи оп­
тимизации не зависит от количества элементов и их геометрических свойств
(считаем, что сетка состоит из однородных элементов), поэтому можем упро­
стить задачу и рассмотреть лишь один элемент в локальной системе координат.
Таким образом, для квадратичного серендипового элемента задача оптимиза­
ции сводится к поиску минимума квадратичной параболы
∫︁1 ∫︁1 ∑︁
8 ∑︁
2
2 2
min 𝑐𝑗 𝑁𝑖,𝑗 (ξ, η)𝑑ξ𝑑η = min 𝐶(27𝑠2 − 12𝑠 + 19) → 𝑠 = , (3.5)
𝑠 𝑠 9
−1 −1 𝑖=0 𝑗=0

где 𝐶 — константа. Исходя из полученной оценки, ожидаемый минимум чис­

ла обусловленности, а также наибольшая скорость сходимости итерационных
методов решения СЛАУ должна быть в окрестности точки 𝑠 = 2/9.
 51

3.5. Основные результаты и выводы по главе 3

1. Представлена общая структура программного комплекса NonLocFEM,

взаимосвязь модулей программы и их предназначение.
2. Предложен способ элементной аппроксимации области нелокального вли­
яния, суть которого заключена в аппроксимации области нелокального
влияния относительно центров элементов. При помощи данного способа
был построен параллельный алгоритм ассемблирования матриц теплопро­
водности и жёсткости.
3. На основе k-d дерева построен быстрый алгоритм аппроксимации области
нелокального влияния. Суть метода заключена в разделении области на
квадратные ячейки, не превышающие радиус нелокальности, для того,
чтобы сузить область поиска и за счёт этого ускорить общее время работы
алгоритма аппроксимации.
4. Рассмотрено семейство базисов квадратичного серендипового элемента;
получена оценка параметра базиса, при которой достигается минимальное
число обусловленности матриц теплопроводности и жёсткости.
 52

Глава 4. Анализ результатов расчётов

4.1. Стратегия исследования и обезразмеривание

В дальнейших расчётах, для изучения качественных различий между

классической (локальной) и нелокальной моделями, проведём процедуру обез­
размеривания основных расчётых параметров уравнений теплопроводности
(1.2) и равновесия (1.5), где безразмерные параметры будем обозначать теми
же символами, но с чертой над ними

𝑥 𝑇 λ
̂︀ 𝐸
𝑥= , 𝑇 = , λ= , 𝐸= , α𝑇 = α𝑇 𝑇0 .
𝐿 𝑇0 λ0 σ0
Здесь 𝐿 — характерный размер области; 𝑇0 — нормализующий множитель для
температуры; λ0 — нормализующий множитель для тензора теплопроводно­
сти; σ0 — нормализующий множитель для напряжений. Также чертой сверху
будем обозначать безразмерные величины, которые образуются посредством
комбинации приведённых выше безразмерных параметров. В расчётах, если не
оговорено иначе, примем безразмерный тензор теплопроводности λ = ̂︀I2 , безраз­
мерный модуль Юнга 𝐸 = 400, коэффициент Пуассона ν = 0.3 и безразмерный
температурный коэффициент линейного расширения α𝑇 = 2.5 · 10−3 .
Стратегия исследования модели подразумевает вариацию основных па­
раметров при фиксации всех остальных. Весовой параметр 𝑝1 и радиус
нелокальности 𝑟 будем варьировать линейно, причём для параметра 𝑝1 будут
рассмотрены четыре сценария: отсутствие нелокальных эффектов (𝑝1 = 1);
локальное слагаемое преобладает над нелокальным (𝑝1 = 0.75); локальное
и нелокальное слагаемые имеют одинаковый вес (𝑝1 = 0.5); и нелокальное
слагаемое преобладает над локальным (𝑝1 = 0.25). Полностью нелокальную
постановку (𝑝1 = 0) рассматривать не будем, так как это приводит к некоррект­
но поставленным краевым задачам, что требует дополнительных рассуждений
при их решении.
 53

Все дальнейшие расчёты будем проводить с использованием квадратич­

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

4.2. Основные особенности решений

Рассмотрим изучение нелокальных моделей теплопроводности

и упругости с решения серии задач на единичном квадрате 𝑆 =
{𝑥 | − 0.5 ⩽ 𝑥1 , 𝑥2 ⩽ 0.5}. Решения будем искать на равномерной сетке 𝑆ℎ
с характерным размером элементов ℎ = 0.004. Поставим граничные и инте­
гральное условия для уравнения теплопроводности
∫︁
𝑛 · 𝑞|𝑥1 =−0.5 = 1, 𝑛 · 𝑞|𝑥1 =0.5 = −1, 𝑇 𝑑𝑆 = 0,
𝑆
а также сформулируем граничные и геометрические условия для уравнения
равновесия

𝑛𝑗 σ𝑗1 |𝑥1 =−0.5 = −1, 𝑛𝑗 σ𝑗1 |𝑥1 =0.5 = 1, 𝑢1 |𝑥1 =0 = 0, 𝑢2 |𝑥2 =0 = 0.

Для определения относительного отклонения будем рассматривать нор­

мированную разность нелокального и локального решений в сечениях вдоль
оси нагружения, где в качестве нормировочного множителя будет выступать
максимальное по модулю значение локального решения
𝑁𝐿 𝐿
𝑇 −𝑇 𝑢𝑁 𝐿
1 − 𝑢𝐿1
𝑇̃︀ = ⃒ ⃒, 𝑢
̃︀1 = ⃒ ⃒.
⃒ 𝐿⃒
max ⃒𝑇 ⃒ max ⃒𝑢𝐿1 ⃒
𝑥∈𝑆 𝑥∈𝑆

Сравнительный анализ для уравнения теплопроводности и равновесия будем

проводить одновременно, так как наблюдаемые в решениях явления весьма
похожи.
 54

Начнём изучение решений с вариации весового параметра модели 𝑝1 . Как

показано на Рис. 4.1, увеличение вклада нелокального влияния увеличива­
ет отклонение решений относительно классического. Причём при параметре
𝑝1 = 0.25 наблюдаем осциляции решения вблизи границ области, где также
наблюдаем кромочный эффект, характеризующийся резкими изменениями ре­
шения.
˜ ˜
T u1
0.03
r = 0.1 r = 0.1
0.04
x1 = 0 x1 = 0 0.02
P 0.02 P
φ = φ2,1 φ = φ2,1 0.01
p1 = 0.75 p1 = 0.75
x1 x1
-0.4 -0.2 0.2 0.4 -0.4 -0.2 0.2 0.4
-0.02 p1 = 0.5 -0.01
p1 = 0.5
p1 = 0.25 -0.02 p1 = 0.25
-0.04
-0.03

а) б)
Рис. 4.1. Решения при вариации весового параметра параметра 𝑝1

Вариация радиуса нелокальности, представленная на Рис. 4.2, также

увеличивает отклонения, но вместе с тем и ширину кромочного эффекта про­
порционально радиусу.
˜ ˜
T u1

p1 = 0.5 0.04 p1 = 0.5

0.02
x1 = 0 x1 = 0
P
φ = φ2,1 0.02 P
φ = φ2,1 0.01
r = 0.05 x1 r = 0.05 x1
-0.4 -0.2 0.2 0.4 -0.4 -0.2 0.2 0.4
-0.02 r = 0.1 -0.01 r = 0.1
r = 0.15
-0.04 r = 0.15 -0.02

а) б)
Рис. 4.2. Решения при вариации радиуса нелокальности 𝑟

Заметим, что в двумерном случае, отклонения также зависят и от рассмат­

риваемого сечения. На Рис. 4.3 представлены распределения температуры и
 55

перемещения в различных сечения и при приближении к свободным от условий

границам отклонения возрастают, так как на них решение также подверже­
но кромочному эффекту. Наибольшие отклонения решений находятся в углах
области, где кромочные эффекты двух границ складываются, и отклонения до­
стигают 0.06 для уравнения теплопроводности и 0.15 для уравнения равновесия
относительно классического решения при параметрах модели 𝑝1 = 0.5 и 𝑟 = 0.1.
˜ ˜
T u1

p1 = 0.5 0.06 p1 = 0.5 0.15

r = 0.1 0.04 r = 0.1 0.10 x 2 = 0.5
P P
φ = φ2,1 0.02 φ = φ2,1 0.05 x 2 = 0.35
x2 = 0
x1 x1
-0.4 -0.2 0.2 0.4 -0.4 -0.2 0.2 0.4
-0.02 -0.05 x2 = 0
x 2 = 0.35
-0.04
-0.10
-0.06 x 2 = 0.5
-0.15

а) б)
Рис. 4.3. Нормированная разница решений в разных сечениях

Помимо уже рассмотренного, решения в нелокальной постановке также

обладают и другими интересными особенностями. Например, в отличие от клас­
сического решения, компонента теплового потока 𝑞 2 не равна нулю вблизи
границ, где задан ненулевой тепловой поток, причём, как показано на Рис. 4.4,
решения обладают симметрией и также зависят от основных параметров моде­
ли. В отличие от температуры, увеличение радиуса нелокальности 𝑟 не влияет
на величину отклонений, но увеличивает размах кромочного эффекта, который
здесь характеризуется увеличением площади областей, окружённых линиями
уровней. Внутри области и на свободных от условий границах компонента плот­
ности теплового потока 𝑞 2 равна нулю.
Для уравнения равновесия в нелокальной постановке ненулевой стано­
вится компонента тензора напряжений σ12 , а заодно и компонента тензора
деформаций ε12 . Её распределение имеет более сложную, но при этом также
 56

q2 q2
x2 x2
0.5 0.5
0.076 0.076

0.25 0.038 0.25 0.038

p1 = 0.5 p1 = 0.5
r = 0.05 r = 0.15
P x1 0 P x1 0
-0.5 -0.25 φ= φ2,1 0.25 0.5 -0.5 -0.25 φ = φ2,1 0.25 0.5

-0.038 -0.038
-0.25 -0.25

-0.076 -0.076
-0.5 -0.5

а) б)
Рис. 4.4. Распределение компоненты 𝑞 2 при вариации 𝑟

симметричную форму и представлено на Рис. 4.5. Здесь, как и в случае с

компонентой плотности теплового потока 𝑞 2 , вариация радиуса практически
не оказывает влияния на максимальные значения напряжения и также увели­
чивает размах линий уровня. Отметим ещё, что на всех границах, включая те,
где заданы нагружения, а также в центре области, значения σ12 равны нулю.
σ 12 σ 12
x2 x2
0.5 0.5
0.0092 0.0096

0.25 0.0046 0.25 0.0048

p1 = 0.5 p1 = 0.5
r = 0.05 r = 0.15
P x1 0 P x1 0
-0.5 -0.25 φ= φ2,1 0.25 0.5 -0.5 -0.25 φ = φ2,1 0.25 0.5

-0.0046 -0.0048
-0.25 -0.25

-0.0092 -0.0096
-0.5 -0.5

а) б)
Рис. 4.5. Распределение компоненты σ12 при вариации 𝑟
 57

4.3. Исследование функций нелокального влияния

Изучим теперь влияние выбора функции нелокального влияния φ. Ранее

уже были определены два параметрических семейства функций: полиномиаль­
ные φ𝑃𝑝,𝑞 (1.9) и экспоненциальные φ𝐸
𝑝,𝑞 (1.11). На их примере покажем как

вариация основных параметров влияет на результаты решений. Не теряя общно­

сти исследования ограничимся решением только уравнения теплопроводности,
но при этом будем сравнивать семейства функций между собой одновременно,
так как их параметры, обозначенные одинаковыми символами, имеют одина­
ковый смысл. Рассматриваемую область и постановку задачи оставим той же,
что уже была рассмотрена в предыдущем разделе.
Начнём исследование с вариации параметра 𝑝. Данный параметр отвечает
за равномерность распределения функции влияния по заданной области 𝑆 ′ (𝑥).
Как показано на Рис. 4.6, отклонения решения растут вместе с параметром 𝑝
и достигают своего максимума, когда функция влияния вырождается в кон­
станту при 𝑝 → ∞. Здесь в качестве области нелокального влияния 𝑆 ′ (𝑥) был
выбран круг с радиусом 0.1, поэтому для экспоненциального семейста функций
был подобран и дисперсионный параметр нелокальности 𝑟 при заданном кван­
тиле 𝑄 = 0.99. Добавим ещё одно замечание, что для того, чтобы поведение
экспоненциального семейства совпадало с полиномиальным, то есть увеличе­
ние параметра 𝑝 приводило к строгому увеличению отклонений, необходимо,
чтобы параметр 𝑝 ⩾ 𝑛.
Увеличение параметра 𝑞 концентрирует распределение функции нелокаль­
ного влияния φ в центре области 𝑆 ′ (𝑥). Соответственно, при его увеличении
отклонения решений относительно классического уменьшаются, что продемон­
стрировано на Рис. 4.7. Здесь важно отметить, что для экспоненциальных
функций дисперсионный параметр нелокальности 𝑟 был выбран на основе
функции с наименьшим значением параметра 𝑞, то есть для всех функций
 58

˜ ˜
T T
0.10
q=1 q = 0.5
n=2 n=2
0.05 0.05

p=1 p = 2, r = 0.033
x1 x1
-0.4 -0.2 0.2 0.4 -0.4 -0.2 0.2 0.4
x 2 = 0.5
p=5 x 2 = 0.5 p = 5, r = 0.070
r = 0.1 -0.05 p=∞ -0.05
p1 = 0.5
p1 = 0.5 p = ∞, r = 0.1

а) б)
Рис. 4.6. Вариация параметра 𝑝

он совпадает. Сделано это по причине того, что данные параметры связаны

между собой и если при изменении параметра 𝑞 подобрать подходящий под
установленный квантиль 𝑄 величину 𝑟, распределения функций для всех 𝑞
будут совпадать.
˜ ˜
T T
p=1 p=2 0.06
0.06 n=2
n=2 0.04
0.04
0.02
0.02 q = 10 q=3
x1 x1
-0.4 -0.2 0.2 0.4 -0.4 -0.2 0.2 0.4
x 2 = 0.5 -0.02 q=3 x 2 = 0.5 -0.02 q=1
r = 0.1 -0.04 q=1 r = 0.033 q = 0.5
-0.04
p1 = 0.5 p1 = 0.5
-0.06 -0.06

а) б)
Рис. 4.7. Вариация параметра 𝑞

Параметр 𝑛 является геометрическим, поэтому его изменение влияет на

форму области 𝑆 ′ (𝑥). Вместе с его увеличением увеличивается и покрываемая
областью нелокального влияния площадь и, как показано на Рис. 4.8, откло­
нение решений. В силу независимости величины дисперсионного параметра 𝑟
от параметра 𝑛 в уравнении (1.12), параметр 𝑟 подбирался по правилу «3 сиг­
ма» на основе длины главной полуоси области 𝑆 ′ (𝑥). Для экспоненциальных
 59

функций при заданных параметрах 𝑝 и 𝑞 влияние параметра 𝑛 становится несу­

щественным, когда он больше 2.
˜ ˜
T T
p=1 p=2 0.06
q=1 0.05 q = 0.5
0.04
0.02
n=1 x1 n=1 x1
-0.4 -0.2 0.2 0.4 -0.4 -0.2 0.2 0.4
x 2 = 0.5 n=2 x 2 = 0.5 -0.02 n=2
r = 0.1 -0.05 n=∞ r = 0.033 -0.04 n=∞
p1 = 0.5 p1 = 0.5 -0.06

а) б)
Рис. 4.8. Вариация параметра 𝑛

Подведём промежуточный итог касательно выбора функции нелокально­

го влияния φ. Было продемонстрировано, что её выбор может повлиять на
степень отклонения решения, но качественных различий между решениями
при различных параметрах функций нет. Также выбор семейства функций не
оказывает существенного влияния на итоговые результаты. В связи с этими об­
стоятельствами, дальнейшие расчёты были проведены только с использованием
квадратичных парабол φ𝑃2,1 при 𝑛 = 2, в силу того, что расчёты с использовани­
ем этой функции требуют наименьшего количества вычислительных ресурсов и
этапы ассемблирования матрицы, вычисления тепловых потоков и напряжений
происходят гораздо быстрее. Конечно, ещё меньше вычислительных операций
требует функция φ ≡ const, но такая функция является предельной и не обла­
дает свойством монотонного убывания из-за чего её использование формально
является некорректным.

4.4. Принципы Сен-Венана и стабильности теплового потока

Изучение свойств нелокальной теплопроводности и упругости следует на­

чинать с простых задач, на примере которых можно исследовать основные
 60

принципы и положения присущие классическим моделям. К таким положениям

можно отнести принцип стабильности теплового потока [27] и принцип Сен­
Венана [58], согласно которым любые интегрально эквивалентные возмущения
поля являются локальными и вызывают одинаковые распределения поля вдали
от точек приложения возмущения. В частности, для уравнения теплопровод­
ности (1.2) таким полем является поле плотности теплового потока 𝑞, а для
уравнения равновесия (1.5) — поле напряжений σ.
̂︀
Для проверки этого высказывания проведём серию расчётов на прямо­
угольной области 𝑆 = {𝑥| − 5 ⩽ 𝑥1 ⩽ 5, −0.5 ⩽ 𝑥2 ⩽ 0.5}, с введённой на ней
равномерной сеткой 𝑆ℎ состоящей из 1500 × 150 элементов. Сформулируем гра­
ничные и интегральное условия для уравнения теплопроводности
∫︁ ∫︁
𝑛 · 𝑞|𝑥1 =−5 = 𝑓 (𝑥2 ), 𝑛 · 𝑞|𝑥1 =5 = −𝑓 (𝑥2 ), 𝑇 𝑑𝑆 = 0,
𝑆

а также сформулируем граничные и геометрические условия для уравнения

равновесия

𝑛𝑗 σ𝑗1 |𝑥1 =−5 = −𝑓 (𝑥2 ), 𝑛𝑗 σ𝑗1 |𝑥1 =5 = 𝑓 (𝑥2 ), 𝑢1 |𝑥1 =0 = 0, 𝑢2 |𝑥2 =0 = 0,

где 𝑓 — функция задающая возмущение поля. Будем рассматривать три инте­

грально эквивалентных варианта возмущения

𝑓1 (𝑥) = 1, 𝑓2 (𝑥) = 2 − 4|𝑥|, 𝑓3 (𝑥) = 4|𝑥|,

приложение которых, в виде теплового потока и давления заданных на левой

и правой границах области 𝑆, изображено на Рис. 4.9.
Вначале рассмотрим распределение компоненты теплового потока 𝑞 1 в ло­
кальной и нелокальной постановках. Для этого на Рис. 4.10 рассмотрим сечения
𝑥2 = 0 и 𝑥2 = 0.5, где можем увидеть, что несмотря на различные тепловые пото­
ки подведённые к левой и правой границам, при удалении от точек приложения
потоки сливаются в единую поверхность, которая в локальном случае является
плоскостью 𝑞 1 = 1, а в нелокальном представляет более сложную поверхность.
 61

n·q=-f1(x 2) x2 n·q=f1(x 2) nj σ j1=-f1(x 2) x2 nj σ j1=f1(x 2)

x1 x1

а) б)
n·q=-f2(x 2) x2 n·q=f2(x 2) nj σ j1=-f2(x 2) x2 nj σ j1=f2(x 2)
x1 x1

в) г)
n·q=-f3(x 2) x2 n·q=f3(x 2) nj σ j1=-f3(x 2) x2 nj σ j1=f3(x 2)
x1 x1

д) е)
Рис. 4.9. Тепловые (а, в, д) и механические (б, г, е) нагружения, прикладывае­
мые к пластине на левой и правой границах

Теперь рассмотрим распределение компоненты тензора напряжений σ11 .

Здесь аналогично ранее рассмотренной компоненте плотности теплового пото­
ка 𝑞 1 решения сливаются в общую поверхность, которая, как и в предыдущем
случае, в локальной постановке является плоскостью σ11 = 1, а в нелокальной
представляет некоторую поверхность. На Рис. 4.11 представлены распределе­
ния компоненты тензора напряжений σ11 в сечениях 𝑥2 = 0 и 𝑥2 = 0.5 в
локальной и нелокальной постановках.
Рис. 4.10 и 4.11 имеют похожие кривые, однако, стабилизация теплового
потока происходит заметно быстрее стабилизации напряжений. Для иллюстра­
ции на Рис. 4.12 представим логарифмическую разницу решений полученных
при нагружениях с использованием функций 𝑓1 и 𝑓2 , где для обобщения
обозначений по оси ординат решения будем обозначать символами 𝒮1 и 𝒮2
соответственно. В полученных распределениях локальные и нелокальные реше­
ния имеют одинаковый характер сходимости, но в центре области расхождения
между ними начинают возрастать, причём в случае тепловых потоков разница
достигает около двух порядков, однако, учитывая порядок величин, такое рас­
хождение можно связать с погрешностями вычислений. Также стоит отметить,
что стабилизация теплового потока происходит монотонно, в то время как ста­
 62

q1 q1
2.0 2.0

x2 = 0 1.5 f = f2 x2 = 0 1.5 f = f2
f = f1 f = f1
p1 = 1 p1 = 0.5
1.0 1.0
r = 0.1
P
0.5 f = f3 φ = φ2,1 0.5 f = f3

x1 x1
-4 -2 2 4 -4 -2 2 4

а) б)
q1 q1
2.0 2.0

x 2 = 0.5 1.5 f = f3 x 2 = 0.5 1.5 f = f3

f = f1
p1 = 1 p1 = 0.5
1.0 1.0 f = f1
r = 0.1
P
0.5 f = f2 φ = φ2,1 0.5
f = f2

x1 x1
-4 -2 2 4 -4 -2 2 4

в) г)
Рис. 4.10. Распределение компоненты плотности теплового потока 𝑞 1 в сечениях
(а, б) 𝑥2 = 0 и (в, г) 𝑥2 = 0.5 в (а, в) локальном и (б, г) нелокальном случаях

билизация напряжений имеет осцилирующий характер. Это можно понять по

изломам графика, то есть в этих точках происходит пересечение кривых.
Теперь рассмотрим сечение вдоль оси O𝑥2 . В этом сечении решения
обладают кромочным эффектом, который проявляется в снижении уровня рас­
сматриваемой величины на свободных границах области и компенсирующим
это снижение повышении этой величины в центре. Вариация радиуса нелокаль­
ного влияния 𝑟 увеличивает размах кромочного эффекта, а вариация весового
параметра 𝑝1 влияет на величину отклонения. При этом заметим, что при фик­
сированном радиусе нелокальности 𝑟 все решения, при различных значениях
параметра 𝑝1 , пересекаются в общих точках. В силу того, что графики компо­
ненты теплового потока 𝑞 1 и компоненты напряжений σ11 в этом сечении не
 63

σ 11 σ 11
2.0 2.0
f = f2 f = f2
x2 = 0 1.5 x2 = 0 1.5
p1 = 1 f = f1 p1 = 0.5
f = f1
1.0 1.0
r = 0.1
P
0.5 φ = φ2,1 0.5
f = f3 f = f3
x1 x1
-4 -2 2 4 -4 -2 2 4

а) б)
σ 11 σ 11
2.0 2.0
f = f3
x 2 = 0.5 x 2 = 0.5 f = f3
1.5 1.5
p1 = 1 f = f1 p1 = 0.5
1.0 1.0
r = 0.1 f = f1
P
0.5 φ = φ2,1 0.5
f = f2 f = f2
x1 x1
-4 -2 2 4 -4 -2 2 4

в) г)
Рис. 4.11. Распределение компоненты тензора напряжений σ11 в сечениях (а, б)
𝑥2 = 0 и (в, г) 𝑥2 = 0.5 в (а, в) локальном и (б, г) нелокальном случаях

отличаются, изобразим их на общем Рис. 4.13, указав по оси ординат обобща­

ющий их символ 𝒮.
Во всех сечениях равнодействующие компоненты теплового потока 𝑞 1 и
напряжения σ11 сохраняются и равны приложенным нагружениям:
∫︁0.5 ∫︁0.5 ∫︁0.5 ∫︁0.5
𝑞 1 𝑑𝑥2 = 𝑓𝑖 (𝑥2 )𝑑𝑥2 , σ11 𝑑𝑥2 = 𝑓𝑖 (𝑥2 )𝑑𝑥2 , 𝑖 = 1,3.
−0.5 −0.5 −0.5 −0.5

Это свидетельствует о выполнении принципов стабильности теплового потока

и Сен-Венана, а также сохранении балансных соотношений.
 64

log10(|1-2|) (x 2 = 0) log10(|1-2|) (x 2 = 0.5)

x1 x1
-4 -2 2 4 -4 -2 2 4
-2

-5 σ NL -4
11
σ L11 σ NL
-6 11

-8
σ L11 -10 q NL
-10 1
q NL
1
q L1 -12
-15 q L1 -14

а) б)
Рис. 4.12. Логарифмическая разность решений при различных нагружениях в
сечениях (a) 𝑥2 = 0 и (б) 𝑥2 = 0.5

 

1.00 1.0
r = 0 r = 0.05 r = 0.1 r = 0.15 p1 = 1 p1 = 0.75 p1 = 0.5 p1 = 0.25
0.95 0.9
p = 0.5 r = 0.1
0.90
x1 = 0 x1 = 0 0.8
P 0.85 P
φ = φ2,1 φ= φ2,1
0.80 0.7
x2 x2
-0.4 -0.2 0.2 0.4 -0.4 -0.2 0.2 0.4

а) б)
Рис. 4.13. Распределение компоненты теплового потока 𝑞 1 и компоненты тензора
напряжений σ11 в сечении 𝑥1 = 0 при вариации (а) 𝑟 и (б) 𝑝1

4.5. Растяжение пластины со ступенчатым переходом

Большой интерес представляет поведение модели на областях с сингуляр­

ными точками, где решения стремятся к бесконечности при дроблении сетки.
К таким областям относятся области со ступенчатыми переходами, часто воз­
никающими в различных конструкциях и деталях. Для изучения особенностей
будет достаточно одного перехода, поэтому рассмотрим простейшую Т-образ­
ную область 𝑆 заключённую в единичный квадрат 𝑆 = {𝑥 | 0 ⩽ 𝑥1 , 𝑥2 ⩽ 1}.
 65

Приложим следующие граничные и геометрические условия

𝑛𝑗 σ𝑗2 |𝑥2 =0 = −1, 𝑢2 |𝑥2 =1 = 0, 𝑢1 |𝑥1 =0.5 = 0.

Область с происллюстрированными граничными условиями и интересующими

сечениями AB, CD, EF и GH представлена на Рис. 4.14. Решение будем искать
на равномерной сетке 𝑆ℎ с характерным размером элементов ℎ = 1/300.

G H

E F
C D

x2 A B
x1

σ j2nj = -1
Рис. 4.14. Т-образная область с приложенной нагрузкой

Распределение полей компоненты деформации ε22 в локальном и нелокаль­

ном случаях при различных радиусах нелокального влияния 𝑟 представлены на
Рис. 4.15. В нелокальном случае линии уровня изменяют свой характер вблизи
границ области, особенно в области ступенчатого перехода, где помимо прочих
искажений наблюдаем увеличение деформаций, причём как положительных,
так и отрицательных. Линии уровня в областях с отрицательной деформацией
становятся более выраженными и их размах увеличивается вместе с радиусом
нелокальности 𝑟. В отличие от классического случая, в нелокальном в точ­
ках приложения нагружения поле деформации ε22 неравномерно, несмотря на
 66

равномерный характер нагружения. Это характеризуется повышенными зна­

чениями деформации в углах кромки, к которой приложена нагрузка. Вблизи
точек закрепления линии уровня терпят излом, характеризующийся резкой сме­
ной направления линии, которая по итогу, в отличие от классического случая,
выходит на границу области не под прямым углом. Крутизна излома умень­
шается при увеличении радиуса нелокальности, а расстояние от границы до
излома увеличивается.
ε22 ε22
x2 x2
1.0 0.0100 1.0
0.0125
0.8 0.8
0.0075 0.0100

0.6 0.6 0.0075

0.0050
0.0050
0.4 p1 = 1 0.4 p1 = 0.5
r =0 0.0025 r = 0.05 0.0025
0.2 P
φ = φ2,1 0.2 P
φ = φ2,1 0
0
x1 x1 -0.0025
0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

а) б)
ε22 ε22
x2 x2
1.0 1.0
0.0125 0.0115
0.8 0.8
0.0100

0.6 0.0075 0.6 0.0069

0.0050
0.4 p1 = 0.5 0.4 p1 = 0.5
0.0023
r = 0.1 0.0025 r = 0.15
0.2 P 0.2 P
φ = φ2,1 φ = φ2,1
0
-0.0023
x1 -0.0025 x1
0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

в) г)
Рис. 4.15. Распределение компоненты тензора деформации ε22 при различных
радиусах нелокальности 𝑟
 67

Аналогичные распределения деформации в областях со ступенчатыми пе­

реходами были получены в результате серии экспериментов под руководством
А.В. Андреева с точным измерением деформаций на пластинах из оптиче­
ски активных материалов [20]. Деформации измерялись в различных сечениях
пластины с использованием 120 тензородатчиков с одинковым механическим
сопротивлением. Также был использован альтернативный метод измерения де­
формации с использованием оптических приборов, где при помощи нанесённой
на пластину сетки измерялись перемещения её узлов. При помощи закона Гука
и экспериментальных данных о деформациях были определены напряжения и
найдены их равнодействующие в сечениях. В работе А.В. Андреева было показа­
но, что на участке, отстоящем от ступенчатого перехода примерно на четверть
ширины ступени, наблюдались наиболее резкие изменения полей деформации
и напряжений. Кроме того, в результате экспериментов было установлено, что
равнодействующая напряжений, расчитанная по формуле
∫︁𝑙2
σ*22 = σ22 𝑑𝑥1 ,
𝑙1

в сечениях, достаточно удалённых от ступенчатого перехода, равна при­

ложенной нагрузке, а результирующие напряжения в сечениях близких к
ступенчатому переходу, оказались в среднем на 30% меньше приложенного на­
гружения. Эффект сохранялся при различных значениях нагрузки, а также
при двух предельных формах закона Гука, соответствующих плоскому дефор­
мированному и плоскому напряжённому состояниям. Объяснением подобного
эффекта в работе А.В. Андреева вероятно является тот факт, что классическая
теория упругости при расчёте напряжений не позволяет учесть структурных
особенностей материала, что и привело к неправильной интерпретации полу­
ченных данных.
Действительно, если подставить значения деформации, представленные
на Рис. 4.15, в классический закон Гука, то получим результаты похожие на
 68

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

пряжения σ*22 , представленные на Рис. 4.16, также демонстрируют снижение
напряжения сразу после ступенчатого перехода, но не на 30%, как это было
описано в экспериментах. В нижней части области ситуация обратная, резуль­
тирующее напряжение выше приложенной нагрузки, а в области стыковки двух
частей области наблюдаем резкий перепад напряжений. Максимум отклонения
находится на нижней и верхней границах области. Отметим, что наибольший
вклад на величину отклонения влизи границ области оказывает весовой па­
раметр 𝑝1 , а на величину отклонения внутри области, а также на ширину
отклонения, радиус нелокального влияния 𝑟. При расчёте напряжений по фор­
муле (1.6) равнодействующая напряжений во всех сечениях равна постоянной,
интегрально совпадающей с приложенной нагрузкой.
σ *22 σ *22
0.70 1.0
x 1 = 0.5 0.9 x 1 = 0.5
0.65
r = 0.1 0.8 r = 0.1
0.60 P
r = 0.15 r = 0.1 φ= φ2,1 0.7 P
φ = φ2,1
0.55 p1 = 0.25
0.6 p1 = 0.5
0.50
r =0 0.5
r = 0.05 p1 = 1 p1 = 0.75
0.45 x2 x2
0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

а) б)
*
Рис. 4.16. Равнодействующая напряжения σ22 при вариации (а) 𝑟 и (б) 𝑝1

На Рис. 4.17 представлены распределения компоненты напряжения σ22

в сечениях, указанных на Рис. 4.14, при различных весовых параметрах 𝑝1
и фиксированном радиусе нелокального влияния 𝑟 = 0.1. Во всех сечениях
увеличение вклада нелокального влияния увеличивает отклонение решения
относительно классического. При этом в нижней части области отклонения
гораздо более заметные, чем в верхней. Как и в задаче о растяжении пря­
моугольной пластины, на примере которой ранее был рассмотрен принцип
 69

Сен-Венана, здесь наблюдаем существенное снижение напряжения σ22 на сво­

бодных от условий границах области. Внутри области напротив происходит
повышение уровня напряжений. Отдельно стоит отметить сечение AB, в ко­
тором у нелокальных решений образуются «горбы» отстоящие на радиус
нелокальсти 𝑟 от границ. В сечениях близких к концентратору CD и EF наблю­
даем существенное снижение максимального уровня напряжений. А в сечении
EF различия в решениях не настолько существенные, как в других сечениях,
но также подчиняются всем ранее описанным особенностям. Дополнительно от­
метим, что при фиксированном радиусе нелокального влияния 𝑟 все решения
пересекаются в общих точках.

σ 22 σ 22
1.3 CD (x 2 = 0.45)
AB (x 2 = 0.1)
1.05
p1 = 0.75 1.2
P
1.00 r = 0.1 φ = φ2,1
p1 = 1 1.1
p1 = 0.5
0.95 1.0
P
r = 0.1 φ = φ2,1 p1 = 0.25
0.90 0.9 p1 = 1
p1 = 0.75
0.85 0.8 p1 = 0.5
p1 = 0.25
x1 x1
0.3 0.4 0.5 0.6 0.7 0.3 0.4 0.5 0.6 0.7

а) б)
σ 22 σ 22
1.5 EF (x 2 = 0.52) 0.8
p1 = 1
p1 = 0.25
GH (x 2 = 0.75)
0.6 p1 = 1
1.0

p1 = 0.5 0.4
p1 = 0.75 r = 0.1 P
φ = φ2,1
0.5 p1 = 0.25
P
r = 0.1 φ = φ2,1 0.2
p1 = 0.75 p1 = 0.5
x1 x1
0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

в) г)
Рис. 4.17. Распределение напряжений σ22 при различных весовых параметрах
𝑝1
 70

4.6. Задача Кирша с обобщением на эллиптические вырезы

Продолжим изучение модели на примере решения задачи Кирша с обобще­

нием на эллиптические вырезы. Для этого рассмотрим область 𝑆, заключённую
в квадратную область 𝑆 = {𝑥 | − 1 ⩽ 𝑥1 , 𝑥2 ⩽ 1}, с эллиптическим вырезом по
центру, где главные оси выреза сонаправлены с осями координат и имеют длины
𝑅1 и 𝑅2 соответственно. Поставим граничные и геометрические условия

𝑛𝑗 σ𝑗1 |𝑥1 =−1 = −1, 𝑛𝑗 σ𝑗1 |𝑥1 =1 = 1, 𝑢1 |𝑥1 =0 = 0, 𝑢2 |𝑥2 =0 = 0.

Схематичное представление области и приложенных нагружений изображены

на Рис. 4.18, где также введена дополнительная угловая координата θ необходи­
мая для анализа распределения интересуемых величин на кромке AB. Решение
будем искать на структурированной сетке 𝑆ℎ , построенной при помощи про­
граммного комплекса Abaqus [1], пример которой при 𝑅1 = 0.2 и 𝑅2 = 0.4 с
характерным размером элементов ℎ = 0.05 представлен на Рис. 4.19. Расчёты
описанные далее были проведены на более подробных сетках с характерным
размером элементов ℎ = 0.005 при различных параметра 𝑅1 и 𝑅2 не превы­
шающих 0.1.
Наибольший интерес исследования представляет распределение дефор­
мации и напряжений на кромке эллиптического выреза, но так как задача
обладает симметрией, сократим рассматриваемую область до дуги AB. Для
удобства исследования введём угловой параметр θ, относительно которого па­
раметризуем координаты дуги следующим образом
⎧
⎨𝑥1 (θ) = 𝑅1 cos θ,
⎪

⎩𝑥2 (θ) = 𝑅2 sin θ,

⎪
 71

B
x2
R2
nj σ j1 = -1

nj σ j1 = 1
θ x1 A
R1

Рис. 4.18. Постановка задачи Кирша

где угол θ принимает значения от 0 до π/2. Затем вычислим длину дуги 𝑙 с

зависимостью от угла θ

∫︁θ
√︃(︂ )︂2 (︂ )︂2
𝜕𝑥1 𝜕𝑥2
𝑙(θ) = + 𝑑φ
𝜕φ 𝜕φ
0

и наконец обезразмерим этот параметр

𝑙(θ)
𝑙(θ) = (︁ π )︁ . (4.1)
𝑙
2
Дальнейшие результаты на кромке AB будем рассматривать в координатах без­
размерного параметра длины 𝑙. В расчётах будем варьировать отношение длин
полуосей выреза таким образом, чтобы максимальная длина была равна 0.1,
т.е. max(𝑅1 , 𝑅2 ) = 0.1. Также для дальнейшего анализа введём величину отно­
шения длин полуосей ρ = 𝑅2 /𝑅1 .
Для данной постановки задачи известно, что максимальные значения
компоненты тензора напряжений σmax
11 находятся в верхней и нижней точках

эллипса и образуют линейную зависимость относительно параметра ρ и при­

 72

x2
x1

Рис. 4.19. Пример конечно-элементной сетки на области с эллиптическим выре­

зом

кладываемого нагружения σ0 [23, 25]

σmax
11 = (1 + 2ρ) σ0 .

Однако в нелокальном случае максимальный уровень напряжения снижается

и начинает зависеть ещё и от весового параметра 𝑝1 . Рассмотрим результаты
представленные в Таб. 4.1, где выписаны значения σmax
11 при различных значе­

ниях параметра ρ и весового параметра 𝑝1 . Заметим, что данные в локальном

случае (𝑝1 = 1) хорошо согласуются с представленной выше зависимостью, а
в нелокальном необходимо добавить дополнительный множитель κ, зависящий
от весового параметра 𝑝1 ,

σmax
11 = κ (1 + 2ρ) σ0 .

Такая зависимость не имеет строгого теоретического доказательства и получе­

на эвристически, однако, она может быть полезна для оценки максимальных
 73

значений напряжений в практических расчётах. Для оценки параметра κ при

фиксированном значеним 𝑝1 необходимо провести 𝑛 расчётов при различных
соотношениях длин полуосей и выполнить осреднение согласно следующей фор­
муле
𝑛
1 ∑︁ σmax
11 (ρ𝑖 , 𝑝1 )
κ(𝑝1 ) = , (4.2)
𝑛 𝑖=0 σmax
11 (ρ𝑖 , 1)

где ρ𝑖 — отношение длин полуосей в 𝑖-ом расчёте. По результатам, представлен­

ным в Таблице 4.1, значения κ линейно зависят от весового параметра 𝑝1 .

Таблица 4.1. Максимальный уровень напряжения σ11 при

вариации отношения длин полуосей ρ и весового
параметра 𝑝1 , где 𝑟 = 0.05
Отношение длин Весовые параметры
полуосей ρ 𝑝1 = 1 𝑝1 = 0.75 𝑝1 = 0.5 𝑝1 = 0.25
0.5 2.012 1.783 1.537 1.510
0.75 2.578 2.235 1.919 1.727
1 3.053 2.696 2.308 1.937
1.25 3.532 3.123 2.670 2.139
1.5 4.012 3.551 3.031 2.404
κ 1 0.881 0.755 0.652

Согласно результатам из Таблицы 4.1, при 𝑝1 = 0.75 и 𝑝1 = 0.5 осредне­

ние параметра κ (4.2) не нужно, так как результаты отношений напряжений в
локальном и нелокальном случаях близки при любых значениях ρ, представ­
ленных в таблице, но при 𝑝1 = 0.25 такая зависимость нарушается, однако,
при осреднении зависимость параметра κ от весового параметра 𝑝1 становится
близкой к линейной. Разумеется все эти рассуждения требуют более деталь­
ного теоретического рассмотрения, так как имеющихся эмпирических данных
недостаточно для утверждения явных зависимостей.
Вместе с напряжением σ11 при увеличении значения ρ увеличивается и
максимальный уровень деформации ε11 . Причём если обратить внимание на
 74

Рис. 4.20 то можем заметить, что помимо повышения максимального уровня

деформации она также становится более сконцентрированной в верхней и соот­
ветственно нижней точках выреза. В нелокальном случае наблюдаем похожий
эффект, однако, здесь также появляется область с отрицательными значения­
ми деформации, которая находится рядом с концентратором. Такой же эффект
можно наблюдать при решении нелокальных задач на других областях с кон­
центраторами [103] и экспериментах [20].
ε11
ρ = 1.5 ε11
0.010 ρ = 1.5
p1 = 1 0.012
0.008 p1 = 0.5
r =0 0.010
r = 0.05
0.006 P 0.008
φ = φ2,1 ρ = 0.75 P
0.006 φ = φ2,1 ρ = 0.75
0.004 ρ = 0.5
ρ = 1.25 0.004 ρ = 0.5
0.002 0.002 ρ = 1.25
ρ=1 ρ=1
l l
0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

а) б)
Рис. 4.20. Распределение компоненты деформации ε11 на кромке AB в (а)
локальном и (б) нелокальном случаях при различных соотношениях длин по­
луосей выреза

При уменьшении параметра 𝑝1 , в отличие от напряжения σ11 , деформация

ε11 увеличивается в зоне концентрации. Также увеличивается и смежная с ней
зона с отрицательными значениями деформации, а сами значения в ней увели­
чиваются по модулю. Увеличение радиуса нелокального влияния 𝑟 оказывает
похожий, но менее выраженный эффект. Результаты представлены на Рис. 4.21.
В дополнение отметим, что вариация 𝑟 практически не оказывает влияния на
величину максимального значения напряжения σ11 .

4.7. Тепловые деформации в областях с эллиптическими вырезами

Рассмотрим задачу на той же области, что и задачу Кирша, но сменим

тип нагружений с механических на тепловые, то есть поставим следующие гра­
 75

ε11 ε11

0.008 ρ=1 0.008 ρ=1

r = 0.1
0.006 r = 0.05 p1 = 1 0.006 p1 = 0.5
r =0
P P
φ = φ2,1 φ = φ2,1
0.004 0.004
0.002 p1 = 0.25
p1 = 0.75 0.002
p1 = 0.5 r = 0.025
l r = 0.05
0.2 0.4 0.6 0.8 1.0 l
0.2 0.4 0.6 0.8 1.0

а) б)
Рис. 4.21. Распределение компоненты деформации ε11 на кромке AB при вари­
ации (а) весовопого параметра 𝑝1 и (б) радиуса нелокальности 𝑟

ничные и геометрические условия

𝑛 · 𝑞|𝑥1 =−1 = 1, 𝑛 · 𝑞|𝑥1 =1 = −1, 𝑢2 |𝑥1 =0 = 0.

Графическое изображение постановки задачи представлено на Рис. 4.22. Также,

для достижения единственности решения, добавим интегральные условия на
искомую температуру и первую компоненту вектора перемещения
∫︁ ∫︁ ∫︁ ∫︁
𝑇 𝑑𝑆 = 0, 𝑢1 𝑑𝑆 = 0.
𝑆 𝑆

Такая постановка удобна тем, что позволяет качественно изучить поведение

температурных напряжений без появления дополнительных напряжений со
стороны возможных концентраторов, обусловленных граничными или геомет­
рическими условиями.
Перед изучением полей температурных напряжений обратим внимание,
что аналогично компоненте тензора напряжения σ11 максимальные значения
компоненты вектора плотности теплового потока 𝑞 max
1 находятся на верхней и
нижней точках эллиптического выреза и они подчинены следующей законо­
мерности

𝑞 max
1 = (1 + ρ)𝑞𝑜 ,
 76

B
x2
R2

n·q = -1
n·q = 1 θ x1 A
R1

Рис. 4.22. Постановка задачи с тепловыми нагружениями на области с эллип­

тическим вырезом

где 𝑞0 — величина подаваемого теплового потока. В нелокальном случае

аналогично напряжениям величина 𝑞 max
1 снижается при увеличении вклада
нелокального влияния, однако, по результатам представленным в Таблице 4.2
не удаётся также легко определить получившуюся зависимость от весового па­
раметра 𝑝1 , так как при весах 𝑝1 = 0.5 и 𝑝1 = 0.25 и значении ρ ⩽ 1 величины
𝑞 max
1 достаточно близки и линейная зависимость от 𝑝1 наблюдается только при
ρ = 1.5.

Таблица 4.2. Максимальное значение компоненты

теплового потока 𝑞 1 при вариации отношения длин
полуосей ρ и весового параметра 𝑝1 , где 𝑟 = 0.05
Отношение длин Весовые параметры
полуосей ρ 𝑝1 = 1 𝑝1 = 0.75 𝑝1 = 0.5 𝑝1 = 0.25
0.5 1.501 1.333 1.316 1.326
0.75 1.753 1.556 1.457 1.462
1 2.001 1.781 1.597 1.591
1.25 2.252 2.002 1.734 1.644
1.5 2.494 2.222 1.927 1.681
 77

Рассмотрим теперь температурные напряжения. Благодаря интегральным

условиям, решения симметричные и все напряжения сконцентрированы вокруг
выреза. Относительно верхней и нижней половин области функции решений
чётные, а относительно левой и правой половин нечётные. Зная эти особенно­
сти, ограничимся изучением распределения полей напряжения лишь на дуге
AB. Для удобства изучения также воспользуемся обезразмеренным парамет­
ром длины 𝑙 (4.1).
При вариации параметра ρ кривые σ11 и σ22 существенно изменяют свои
формы. На Рис. 4.23 представлены распределения этих кривых при различных
параметрах ρ в классическом случае (𝑝1 = 1). У кривой σ11 вместе с увеличени­
ем ρ пиковое значение смещается право вдоль оси O𝑙. Помимо этого оно также
растёт в абсолютных значениях. Максимальная величина σ22 также меняется,
но основная закономерность заключается в увеличении площади под кривой
при росте ρ. Важно отметить, что для обеих кривых в точке 𝑙 = 1 их величины
равны 0, так как в ней происходит смена знака.
σ 11 σ 22
0.05 ρ = 1.5
p1 = 1 0.10
p1 = 1
0.04
r =0 0.08 r =0
ρ = 1.25
0.03 0.06
0.02 ρ=1 ρ = 1.5
0.04
ρ = 0.75
ρ = 0.75 ρ = 1.25
0.01 0.02 ρ = 0.5 ρ=1
ρ = 0.5
l l
0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

а) б)
Рис. 4.23. Распределение напряжения (a) σ11 и (б) σ22 на кромке AB при вари­
ации соотношения длин полуосей

Другими словами, сужение эллиптического выреза вдоль оси O𝑥1 приво­

дит к смещению пиковых значений компоненты напряжений σ11 к верхней и
нижней точкам выреза, а также увеличению их абсолютных значений. Пиковое
 78

значение компоненты напряжений σ22 напротив уменьшается, но при этом зани­

маемые площади с ненулевыми напряжениями увеличиваются. Таким образом,
эллиптические вырезы большая ось которых располагается вдоль линии тока
теплового потока дают меньшие напряжения, чем тем, у которых большая ось
располагается поперёк.
Учёт нелокальных свойств среды приводит к снижению напряжений σ11
и σ22 . Как и во всех расчётах до этого, вариация параметра 𝑟 влияет лишь
на форму распределения полей напряжений, а вариация весового параметра 𝑝1
на величину отклонения. Распределение полей σ11 и σ22 вдоль кривой AB при
вариации 𝑝1 представлено на Рис. 4.24.
σ 11 σ 22
0.04 p1 = 0.75 0.10
ρ=1 p1 = 1 ρ=1
0.08 p1 = 1
0.03 r = 0.05 r = 0.05
φ= P
φ2,1 0.06 P
φ = φ2,1
p1 = 0.75
0.02
p1 = 0.5
0.04
p1 = 0.25
0.01 p1 = 0.25 p1 = 0.5
0.02

l l
0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

а) б)
Рис. 4.24. Распределение напряжения (a) σ11 и (б) σ22 на кромке AB при вари­
ации весового параметра 𝑝1

4.8. Основные результаты и выводы по главе 4

1. Исследованы основные особенности получаемых распределений темпера­

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

2. Исследованы полиномиальные и экспоненциальные семейства функций

нелокального влияния. Показано, что вариация основных параметров
функций не даёт качественных различий между решениями, в связи с
чем для остальных расчётов было решено выбрать в качестве функции
нелокального влияния квадратичную параболу.
3. Продемонстрирована применимость принципов Сен-Венана и стабильно­
сти теплового потока в контексте нелокальных постановок уравнений
равновесия и стационарной теплопроводности соответственно. Показано,
что вдали от точек приложения нагружений, поля напряжений и плот­
ности теплового потока сливаются в единую поверхность, независимо от
вида нагружения. Эти поверхности обладают кромочным эффектом на
свободных от нагружений границах, который характеризуется снижени­
ем значения рассматриваемой величины.
4. На примере растяжения Т-образной пластины проиллюстрировано сни­
жение роли концентратора напряжений в нелокальной постановке, по
сравнению с классической.
5. Исследована задача Кирша с обобщением на эллиптические вырезы. В
нелокальной постановке установлена связь между максимальным уровнем
напряжения, приложенного нагружения и формы выреза.
6. Исследована задача о температурных напряжениях на области с эллип­
тическим вырезом. Результаты демонстрируют существенное снижение
уровня напряжения в концентраторах.
 80

Глава 5. Анализ эффективности программного комплекса

NonLocFEM

5.1. Тестирование алгоритма ассемблирования матриц

Оценим масштабируемость полученного алгоритма ассемблирования мат­

риц теплопроводности и жёсткости (3.3). Для этого проведём серию расчётов
на вычислительном кластере, в состав которого входит 6 вычислительных уз­
лов. На каждом узле кластера установлен 18-ядерный процессор Intel Core i9
10980XE и 128 ГБ оперативной памяти DDR4. Будем рассматривать тестовые
задачи на прямоугольной области 𝑆 = {𝑥| − 5 ⩽ 𝑥1 ⩽ 5, −0.5 ⩽ 𝑥2 ⩽ 0.5}, с
введённой на ней равномерной сеткой 𝑆ℎ . Как и в предыдущей главе, все вы­
числения проведены с использованием квадратичных серендиповых элементов,
интегрирование выполнено гауссовыми квадратурами 3-го порядка (9 квадра­
турных узлов), а в качестве функции нелокального влияния выбрана φ = φ𝑃2,1 ,
выбор которой также обоснован в предыдущей главе. Хранение разреженных
матриц организовано в формате CSR [47], где для хранения индексов ис­
пользовались 64-битные целые числа, а для хранения коэффициентов числа с
плавающей точкой двойной точности. Также была учтена симметрия матрицы,
то есть для хранения и вычисления использовалась только верхняя половина
матрицы.
Начнём исследование масштабируемости алгоритмов на машинах с об­
щей памятью, то есть задействуем только один узел кластера и технологию
параллельного программирования OpenMP. Для этого проведём серию расчё­
тов варьируя разбиение сетки ℎ и радиус нелокальности 𝑟, который в данном
исследовании совпадает с радиусом поиска. Результаты, представленные в Таб­
лице 5.1, свидетельствуют о том, что для хранения матрицы жёсткости K
̂︀ 𝐸

требуется в 4 раза больше оперативной памяти, чем для матрицы теплопровод­

ности K
̂︀ 𝑇 на аналогичной сетке, что было ожидаемо, учитывая размерности
 81

блоков матриц теплопроводности (2.9) и жёсткости (2.10). В обоих случаях,

темпы роста занимаемой оперативной памяти относительно разбиения сетки
линейные в классическом случае и квадратичные в нелокальном. При сравне­
нии времён сборки матриц можем отметить, что в классическом случае время
ассемблирования матрицы жёсткости K
̂︀ 𝐸 приблизительно в 3 раза дольше, чем

матрицы теплопроводности K
̂︀ 𝑇 , однако, в нелокальном случае эта разница со­

кращается до 1.2, чего удалось достичь за счёт оптимизации вычислений с

использованием блочного подхода в ассемблировании.

Таблица 5.1. Требуемая оперативная память и затрачиваемое время при

ассемблировании матриц теплопроводности и жёсткости с
использованием технологии OpenMP
Количество Радиус Требуемая Время, с Время, с
элементов поиска оперативная память 1 поток 18 потоков
K
̂︀ 𝑇 K̂︀ 𝐸 K̂︀ 𝑇 K̂︀ 𝐸 K̂︀ 𝑇 K̂︀ 𝐸
400 × 40 0 6.2 Мб 23.8 Мб 0.061 0.187 0.028 0.048
800 × 80 0 24.5 Мб 95 Мб 0.521 1.701 0.116 0.274
1600 × 160 0 97.8 Мб 379 Мб 5.17 16.76 0.705 2.366
400 × 40 0.05 61.4 Мб 244 Мб 14.4 17.06 1.055 1.288
800 × 80 0.05 548 Мб 2.2 Гб 157 187 11.37 13.3
1600 × 160 0.05 5.5 Гб 22 Гб 1825 2166 132.9 155
400 × 40 0.1 133 Мб 532 Мб 38 44.9 2.68 3.26
800 × 80 0.1 1.34 Мб 5.4 Гб 450 529 32 37.8
1600 × 160 0.1 17 Гб 68 Гб 6134 7266 447 521

Обращаясь к результатам из Таблицы 5.1, можем построить диаграммы

эффективности распараллеливания алгоритмов ассемблирования матриц тепло­
проводности и жёсткости. Исходя из полученных результатов, представленных
на Рис. 5.1, отметим достаточно хорошую эффективность распараллеливания,
которая в нелокальном случае достигает 14 раз при использовании всех 18
ядер процессора. Однако ускорение в классическом случае не такое высокое,
что можно объяснить обобщённостью используемых алгоритмов формирова­
 82

ния портрета матрицы, которые недостаточно эффективно работают на малых

объёмах данных.
t1/t18 t1/t18
14.2 14.1 14.1 14. 13.8 14. 13.9
14 13.6 13.8 13.7 13.7 14 13.2

12  12 
KT KE
10 10

1600×160, r = 0.05

1600×160, r = 0.05
8 7.33 8 7.08

1600×160, r = 0.1

1600×160, r = 0.1
400×40, r = 0.05

800×80, r = 0.05

400×40, r = 0.05

800×80, r = 0.05
6.21

400×40, r = 0.1

800×80, r = 0.1

400×40, r = 0.1

800×80, r = 0.1
6 6
4.49
3.9
1600×160

1600×160
4 4
800×80

400×40

800×80
2.18
2 2
400

0 0

а) б)
Рис. 5.1. Эффективность распараллеливания алгоритма сборки матриц (а) теп­
лопроводности и (б) жёсткости при использовании технологии OpenMP

Теперь исследуем масштабируемость алгоритма ассемблирования матриц

на примере использования технологии MPI. Здесь сконцентрируем внимание
на алгоритме балансировки данных между процессами, так как вопрос хране­
ния матриц в классе нелокальных задач стоит наиболее остро. Для этого будем
рассматривать ту же область, но возьмём более подробную сетку 𝑆ℎ , состоя­
ющую из 3200 × 320 элементов. Радиус нелокальности 𝑟 выберем равным 0.1.
С целью исключения дублирования результатов, остановимся на рассмотрении
лишь матрицы теплопроводности.
Проведём сравнение двух запусков на 6 вычислительных узлах кластера.
Первый запуск осуществим без балансировки данных, то есть распределения уз­
лов сетки между процессами оставим равномерным. Второй запуск осуществим
с балансировкой данных между процессами, которая выполняется по алгоритму
описанному в третьей главе. На Рис. 5.2 представлены полученные распреде­
ления затрачиваемой оперативной памяти для хранения блоков матрицы на
каждом узле. В варианте без балансировки данные распределены не равномер­
но и некоторым узлам кластера достался объём данных больше, чем другим. В
 83

варианте с балансировкой данные распределены, в рамках допустимой погреш­

ности, равномерно. Из чего можно сделать вывод, что алгоритм балансировки
данных работает корректно.

36.9 Gb
63.1 Gb

60.7 Gb 36.9 Gb 36.9 Gb

23 Gb
36.9 Gb 36.9 Gb
26.4 Gb
23.1 Gb
25.1 Gb 36.9 Gb

а) б)
Рис. 5.2. Распределение размеров блоков матрицы теплопроводности между 6
процессами в случае (а) без балансировки и (б) с балансировкой данных

В качестве замечания можно ещё отметить, что алгоритм балансировки

данных также ускоряет работу метода сопряжённых градиентов, так как опе­
рация умножения матрицы на вектор становится одинаковой по трудозатратам
на каждом узле кластера. Более подробно про другие возможные методы уско­
рения сходимости решателей СЛАУ обсудим в последующих разделах.

5.2. Анализ скорости сходимости при оптимизации базиса

конечных элементов

В работе повсеместно использовались квадратичные серендиповые элемен­

ты, в связи с чем возникла потребность оптимизиации их базиса для ускорения
сходимости итерационных методов решения СЛАУ. Отметим, что у квадратич­
ного серендипового элемента есть целое параметрическое семейство базисов,
 84

представленное в третьей главе диссертации. В той же главе была получена

оценка оптимального значения параметра 𝑠 (3.5) согласно которой скорость
сходимости должна быть максимальной в точке 𝑠 = 2/9. Проверим эту оценку
на практике, а заодно исследуем влияние других параметров модели.
Для проверки корректности оценки решим серию задач варьируя па­
раметр 𝑠. Расчёты снова будем проводить на прямоугольной области 𝑆 =
{𝑥 | − 5 ⩽ 𝑥1 ⩽ 5, −0.5 ⩽ 𝑥2 ⩽ 0.5}, с введённой на этой области равномерной
сеткой 𝑆ℎ = 1000 × 100. Обезразмерим уравнения (1.2) и (1.5) относительно
коэффициента теплопроводности λ и модуля Юнга 𝐸 соответственно

𝑞 σ
̂︀
𝑞= , σ= ,
λ 𝐸

так как параметры λ и 𝐸 в данном случае выступают в роли масштабирующих

множителей для собственных чисел матриц K
̂︀ 𝑇 и K
̂︀ 𝐸 соответственно. Коэф­

фициент Пуассона ν примем равным 0.3. Поставим граничные условия для

уравнения теплопроводности
∫︁
𝑛 · 𝑞|𝑥1 =−5 = −1, 𝑛 · 𝑞|𝑥1 =5 = 1, 𝑇 𝑑𝑆 = 0,
𝑆

и для уравнения равновесия

𝑛𝑗 σ𝑗1 |𝑥1 =−5 = −1, 𝑛𝑗 σ𝑗1 |𝑥1 =5 = 1, 𝑢1 |𝑥1 =0 = 0, 𝑢2 |𝑥2 =0 = 0.

Решение СЛАУ будем искать методом сопряжённых градиентов без использо­

вания каких-либо дополнительных предобуславливателей.
Для определения числа обусловленности (3.4) необходимо вычислить мак­
сильное и минимальное собственные числа матрицы. Для их вычисления к
программному комплексу NonLocFEM была подключена библиотека Spectra
[15], которая позволяет вычислять собственные числа разреженных матриц и
совместима с библиотекой линейной алгебры Eigen [7].
 85

В ходе экспериментов было установлено, что минимальное собственное

число λmin не зависит от параметра базиса 𝑠 и вариаций функции нелокаль­
ности φ, но при этом имеет зависимость от вклада нелокального влияния и
радиуса нелокальности 𝑟, где при их увеличении величина собственного числа
λmin начинает уменьшаться. Однако исходя из результатов, представленных в
Таблицах 5.2 и 5.3, можем сделать вывод, что зависимость от параметров мо­
дели весьма слабая, так как величина изменений не достигает и десятой доли
от числа полученного в классическом случае, то есть минимальное собственное
число не оказывает серьёзного влияния на число обусловленности. Также мож­
но отметить, что минимальные собственные числа матриц теплопроводности
и жёсткости при одинаковых параметрах модели достаточно близки по зна­
чениям.
Таблица 5.2. Зависимость минимального собственного числа
матрицы теплопроводности K ̂︀ 𝑇 от параметров модели 𝑝1 и 𝑟
𝑝1
1 0.75 0.5 0.25
𝑟
0 3.2637 · 10−6 — — —
0.05 — 3.2429 · 10−6 3.2181 · 10−6 3.1793 · 10−6
0.1 — 3.2231 · 10−6 3.1750 · 10−6 3.0997 · 10−6
0.15 — 3.2032 · 10−6 3.1319 · 10−6 3.0218 · 10−6

Таблица 5.3. Зависимость минимального собственного числа

матрицы жёсткости K ̂︀ 𝐸 от параметров модели 𝑝1 и 𝑟
𝑝1
1 0.75 0.5 0.25
𝑟
0 3.2614 · 10−6 — — —
0.05 — 3.2471 · 10−6 3.2328 · 10−6 3.1958 · 10−6
0.1 — 3.2336 · 10−6 3.1911 · 10−6 3.1191 · 10−6
0.15 — 3.2183 · 10−6 3.1489 · 10−6 3.0435 · 10−6
 86

Максимальные собственные числа λmax матриц теплопроводности K

̂︀ 𝑇 и

жёсткости K
̂︀ 𝐸 напротив имеют достаточно сильную зависимость как от пара­

метра базиса 𝑠, так и от весового параметра модели 𝑝1 . При этом зависимость

от весового параметра 𝑝1 удаётся установить эмпирически, она в точности рав­
на λ𝑁 𝐿 𝐿
max (𝑠) = 𝑝1 λmax (𝑠). Зависимость от радиуса нелокальности 𝑟 и функции

нелокальности φ у максимальных собственных чисел отсутствует. Результаты

представлены в Таблицах 5.4 и 5.5, где также представлены данные о количе­
ствах итераций метода сопряжённых градиентов.

Таблица 5.4. Зависимости количества итераций 𝑁 и максимального

собственного числа λmax матрицы теплопроводности K ̂︀ 𝑇 от вариации
параметров 𝑝1 и 𝑠
Значение 𝑝1 = 1 𝑝1 = 0.75 𝑝1 = 0.5 𝑝1 = 0.25
параметра 𝑠 𝑁 λmax 𝑁 λmax 𝑁 λmax 𝑁 λmax
-1/3 6961 15.997 8724 11.997 5036 7.9980 3549 3.9984
-2/9 5600 11.198 5110 8.3986 4341 5.5987 2948 2.7990
-1/9 4540 7.4660 4095 5.5993 3367 3.7327 2360 1.8663
0 3868 5.3332 3640 3.9997 2891 2.6662 1961 1.3328
1/9 3885 5.3332 3442 3.9997 2861 2.6662 1947 1.3328
2/9 4000 5.3332 3600 3.9997 2826 2.6662 2003 1.3328
1/3 3889 5.3332 3567 3.9997 2778 2.6662 1948 1.3328
4/9 3786 5.3332 3474 3.9997 2835 2.6662 1982 1.3328
5/9 4719 7.4660 4415 5.5993 3525 3.7327 2353 1.8663
2/3 5761 11.198 6887 8.3986 4134 5.5987 2982 2.7990
7/9 7218 15.997 6349 11.997 5285 7.9980 3720 3.9984

Теперь рассмотрим данные из таблиц 5.2-5.5 в графическом представ­

лении. На Рис. 5.3 представлены зависимости числа обсусловленности и
количества итераций метода сопряжённых градиентов от параметра базиса 𝑠
и весового параметра модели 𝑝1 . Обратим внимание, что кривые числа обу­
словленности симметричные относительно точки 𝑠 = 2/9, а также, что число
обусловленности снижается по мере уменьшения весового параметра 𝑝1 . Также
 87

Таблица 5.5. Зависимости количества итераций 𝑁 и максимального

собственного числа λmax матрицы жёсткости K ̂︀ 𝐸 от вариации
параметров 𝑝1 и 𝑠
Значение 𝑝1 = 1 𝑝1 = 0.75 𝑝1 = 0.5 𝑝1 = 0.25
параметра 𝑠 𝑁 λmax 𝑁 λmax 𝑁 λmax 𝑁 λmax
-1/3 7454 13.061 6549 9.7954 4758 6.5298 3893 3.2650
-2/9 6261 9.6366 5066 7.2272 4402 4.8178 3344 2.4090
-1/9 4805 7.0840 4820 5.3128 3992 3.5416 2870 1.7707
0 4797 5.4194 3996 4.0643 3053 2.7094 2509 1.3545
1/9 4380 4.5195 3320 3.3897 2703 2.2602 2289 1.1306
2/9 4272 4.2968 3296 3.2222 2727 2.1476 2229 1.0730
1/3 4270 4.2970 3551 3.2224 2641 2.1477 2229 1.0731
4/9 4085 4.3423 3772 3.2566 3124 2.1709 2242 1.0852
5/9 4446 6.1015 3911 4.5760 3704 3.0504 2659 1.5249
2/3 6138 8.8614 4814 6.6458 4464 4.4302 3203 2.2147
7/9 7076 12.436 5567 9.3273 4634 6.2178 3802 3.1083

стоит отметить, что на интервале 0 ⩽ 𝑠 ⩽ 4/9 кривые числа обусловленно­

сти выходят на плато, центром которого по прежнему является точка 𝑠 = 2/9.
Кривые зависимости числа итераций хорошо коррелируют с кривыми числа
обусловленности и минимум числа итераций находится в окрестностях точки
𝑠 = 2/9, а увеличение вклада нелокального влияния ускоряет сходимость ме­
тода при заданных параметрах.
Результаты, представленные на Рис. 5.4, для матрицы жёсткости анало­
гичны результатам, представленным на Рис. 5.3, для матрицы теплопроводно­
сти, однако, здесь графики зависимости числа обусловленности от параметра
𝑠 уже не являются симметричными и плато на них уже менее выражено, а
графики зависимости числа итераций хуже коррелируют с графиками числа
обусловленности. Тем ни менее, минимум числа обусловленности и числа ите­
раций по прежнему находится в окрестностях точки 𝑠 = 2/9 из чего можно
сделать вывод, что оценка (3.5) пригодна для практического использования.
 88

Также важно отметить, что несмотря на то, что количество итераций

в нелокальном случае становится меньше, время затрачиваемое на решение
СЛАУ в нелокальном случае может быть на несколько порядков больше,
чем для аналогичной задачи в классическом случае. Это связано с объёмами
данных, которые занимают матрицы в нелокальном случае и главным сдер­
живающим фактором в скорости решения СЛАУ выступает скорость работы
оперативной памяти, а не процессора.

Cond(K T ) N
p1 = 1 p1 = 1
2000 p1 = 0.75 8000
p1 = 0.75
p1 = 0.5 p1 = 0.5
6000
1500 p1 = 0.25 p1 = 0.25

4000
1000

2000
s s
1 2 1 1 2 1 4 5 2 7 1 2 1 1 2 1 4 5 2 7
- - - 0 - - - 0
3 9 9 9 9 3 9 9 3 9 3 9 9 9 9 3 9 9 3 9

а) б)
Рис. 5.3. Зависимость (а) числа обусловленности и (б) количества итераций 𝑁
для матрицы теплопроводности K ̂︀ 𝑇 от параметров 𝑠 и 𝑝1 при 𝑟 = 0.1


Cond(K E ) N
p1 = 1 8000 p1 = 1
2000
p1 = 0.75 p1 = 0.75
p1 = 0.5 6000 p1 = 0.5
1500
p1 = 0.25 p1 = 0.25

4000
1000

2000
s s
1 2 1 1 2 1 4 5 2 7 1 2 1 1 2 1 4 5 2 7
- - - 0 - - - 0
3 9 9 9 9 3 9 9 3 9 3 9 9 9 9 3 9 9 3 9

а) б)
Рис. 5.4. Зависимость (а) числа обусловленности и (б) количества итераций 𝑁
для матрицы жёсткости K ̂︀ 𝐸 от параметров 𝑠 и 𝑝1 при 𝑟 = 0.1
 89

5.3. Предобуславливание и выбор начального приближения

Оптимизированный базис квадратичных серендиповых элементов позво­

лил ускорить сходимость метода сопряжённых градиентов при решении СЛАУ
более чем в 1.5 раза по сравнению со случаем, когда используется классиче­
ский базис. Однако существует возможность добиться ещё большей скорости
сходимости за счёт использования предобуславливателей или более подходя­
щих начальных условий.
При разработке предобуславливателя важно учитывать специфику за­
дачи, так как универсальных способов эффективного предобуславливания
не существует. Также важно учитывать возможности современных вычисли­
тельных машин: параллельные вычисления зачастую дают заметно больший
выигрыш во времени, в отличие от предобуславливания в силу того, что мно­
гие популярные методы предобуславливания используют обратный алгоритм
Гаусса, который не всегда возможно эффективно распараллелить, а накладные
расходы при этом могут значительно увеличить цену одной итерации.
Пожалуй основной спецификой рассматриваемого класса уравнений яв­
ляется их матрично-векторное представление (2.7) и (2.8), где матричные
̂︀ 𝐿 и нелокаль­
выражения представлены в виде взвешенных сумм локальных K ℱ
̂︀ 𝑁 𝐿 матриц. При этом локальные слагаемые обалают заметно меньшей
ных K ℱ

плотностью заполнения, чем нелокальные. Учитывая эту специфику, а также

уже рассмотренный анализ границ спектров матриц K
̂︀ ℱ , можем построить пре­

добуславливатель используя для этого лишь локальное слагаемое.

Как правило для построения предобуславливателей необходимы данные
о максимальных и минимальных собственных числах, однако, процесс на­
хождения границ спектра, даже для локальной матрицы, является слишком
медленным. Поэтому примем во внимание полученные знания о связи спектров
локальной и полной матриц и в качестве предобуславливателя возьмём непол­
 90

̂︀ 𝐿 , так как оно не требует

ное разложение Холецкого локальной матрицы K ℱ

непосредственного поиска собственных чисел.

Теперь проведём серию расчётов и на её основе определим эффективность
использования предобуславливателя. В этом же исследовании проверим гипоте­
зу о выборе начального приближения, где в качестве начального приближения
̂︀ 𝐿 · X0 = F, где F — вектор правой части,
выберем результат решения СЛАУ K ℱ

на основе которого будем решать полную СЛАУ и X0 — вектор искомого на­

чального приближения. Для теста будем рассматривать задачи из предыдущего
раздела, а базис элементов выберем оптимальным с параметром 𝑠 = 2/9.
По результатам, представленным в Таблицах 5.6 и 5.7 для уравнений
теплопроводности и равновесия соотвественно, можем сделать вывод, что ис­
пользование предобуславливателя позволяет ускорить решение СЛАУ в N раз.
При этом весовые параметры модели на это не оказывают серьёзного вли­
яния. Вместе с этим можно сделать вывод касательно выбора начального
приближения. Однако предложенная гипотеза не даёт ожидаемого эффекта,
количество итераций сокращается несущественно, а время затрачиваемое на
решение СЛАУ для классической задачи не всегда меньше получаемого вы­
игрыша. В качестве замечания стоит добавить, что для краткости записи
(︁ 𝐿 )︁
формулировка в таблице ILLT K ̂︀
𝑇 означает использование предобуславли­
(︁ 𝐿 )︁
вателя, а ILLT K̂︀
𝑇 + 𝑋0 подразумевает комбинацию предобуславливания и

начального приближения.

Таблица 5.6. Количество итераций и затрачиваемое время при

решении СЛАУ уравнения теплопроводности
(︁ 𝐿 )︁ (︁ 𝐿 )︁
𝑝1 Без предобуславливания ILLT K ̂︀
𝑇 ILLT K̂︀
𝑇 + 𝑋0
𝑁 𝑡, с 𝑁 𝑡, с 𝑁 𝑡, с
0.75 5364 2970 2111 1181 2086 1238 (1167)
0.5 4558 2527 1740 972 1742 1045 (974)
0.25 3020 1673 1250 657 1265 779 (708)
 91

Таблица 5.7. Количество итераций и затрачиваемое время при

решении СЛАУ уравнения равновесия
(︁ 𝐿 )︁ (︁ 𝐿 )︁
𝑝1 Без предобуславливания ILLT K ̂︀
𝑇 ILLT K̂︀
𝑇 + 𝑋0

𝑁 𝑡, с 𝑁 𝑡, с 𝑁 𝑡, с
0.75 6330 13238 2902 6217 2736 6223 (5967)
0.5 5167 10811 2290 4915 2389 5442 (5203)
0.25 3779 7888 1718 3695 1740 4032 (3793)

5.4. Основные результаты и выводы по главе 5

1. Исследована масштабируемость программного комплекса NonLocFEM на

машинах с общей и распределённой памятью. Представленные результаты
демонстрируют хорошую эффективность распараллеливания, выполнен­
ную средствами OpenMP, и балансировку данных между процессами, при
использовании технологии MPI.
2. Проведён анализ скорости сходимости метода споряжённых градиентов
при вариации базиса квадратичного серендипового элемента. Показано,
что оценка, предложенная в разделе 3.4, даёт корректный результат и
минимум числа обусловленности матриц теплопроводности и жёсткости, а
также наибольшая скорость сходимости метода сопряжённых градиентов,
находятся в окрестности точки 𝑠 = 2/9.
3. Предложенный способ предобуславливания на основе неполного разло­
жения Холецкого локальных матриц, продемонстрировал существенный
прирост в скорости сходимости метода сопряжённых градиентов, при
решении уравнений теплопроводности и равновесия в нелокальных поста­
новках.
 92

Общие выводы и заключение

1. Рассмотрена иерархия моделей нелокальной теплопроводности и термо­

упругости, предложено и проанализировано два семейства возможных
функций нелокального влияния, заданных на областях, ограниченных
кривыми Ламэ.
2. Разработан численный алгоритм решения интегро-дифференциальных
уравнений на основе метода конечных элементов, проведена работа над
его оптимизацией и подготовкой к использованию в параллельной среде
вычислений.
3. Разработан собственный программный комплекс NonLocFEM, в рам­
ках которого реализованы все предложенные алгоритмы; параллельные
реализации алгоритмов задействуют технологии параллельного програм­
мирования OpenMP и MPI, все исследования и расчёты проведены в
рамках программного комплекса.
4. Проведён качественный анализ сравнения классических теорий теплопро­
водности и термоупругости с их нелокальными постановками, полученные
результаты свидетельствуют о снижении роли концентраторов в распреде­
лениях полей напряжений и плотности теплового потока и в возникнове­
нии кромочных эффектов на свободных от граничных условий границах,
а также определены основные зависимости отклонений нелокальных ре­
шений относительно классическим путём вариации параметров модели.
5. Исследован вопрос сходимости итерационных методов решения СЛАУ
применительно к задачам в нелокальных постановках, предложены спо­
собы ускорения сходимости с применением альтернативных базисов ко­
нечных элементов и предобуславливателей.
 93

Список литературы

1. Abaqus. URL: [Link] (дата обр.

19.07.2024).
2. Ansys. URL: [Link] (дата обр. 19.07.2024).
3. C++ Reference. URL: [Link] (дата обр. 19.07.2024).
4. CMake. URL: [Link] (дата обр. 19.07.2024).
5. Conan. URL: [Link] (дата обр. 19.07.2024).
6. [Link]. URL: [Link] (дата обр. 19.07.2024).
7. Eigen. URL: [Link] (дата обр. 19.07.2024).
8. FEniCS. URL: [Link] (дата обр. 19.07.2024).
9. FreeFEM. URL: [Link] (дата обр. 19.07.2024).
10. JSON. URL: [Link] (дата обр. 19.07.2024).
11. Lohmann N. nlohmann/json. URL: [Link] (дата обр.
19.07.2024).
12. MPI. URL: [Link] (дата обр. 19.07.2024).
13. OpenMP. URL: [Link] (дата обр. 19.07.2024).
14. Paraview. URL: [Link] (дата обр. 19.07.2024).
15. Spectra. URL: [Link] (дата обр. 19.07.2024).
16. TFlex. URL: [Link] (дата обр. 19.07.2024).
17. Абгарян К. К. Многомасштабное моделирование в задачах структурного
материаловедения: монография. Москва : МАКС Пресс, 2017. 284 с. ISBN
978-5-317-05707-7.
18. Абрамовиц М., Стиган И. Справочник по специальным функциям с
формулами, графиками и математическими таблицами / под ред. П. с ан­
гл. под ред. В.А. Диткина и Л.Н. Карамзиной. Москва : Наука, 1979. 832 с.
19. Александреску А. Современное проектирование на C++. Москва : Изда­
тельский дом «Вильямс», 2008. 336 с. ISBN 978-0201704310.
 94

20. Андреев А. В. Инженерные методы определения концентрации напряже­

ний в деталях машин. Москва : Машиностроение, 1976. 72 с.
21. Аэро Э. Л., Кувшинский Е. В. Континуальная теория асимметрической
упругости // Физика твердого тела. 1964. Т. 10, № 9. С. 2689—2699.
22. Аэро Э. Л., Кувшинский Е. В. Основные уравнения теории упругости
сред с вращательным взаимодействием частиц // Физика твердого тела.
1960. Т. 2, № 7. С. 1399—1409.
23. Безухов Н. И. Основы теории упругости и пластичности. Москва : Изда­
тельство «Высшая школа», 1968. 512 с.
24. Белов П. А., Лурье С. А. Векторная градиентная теория упругости //
Композиты и наноструктуры. 2023. С. 1—15. DOI: 10.36236/1999- 7590-
2022-14-1-1-15.
25. Биргер И. А., Шорр Б. Ф., Иосилевич Г. Б. Расчет на прочность деталей
машин: Справочник. 4-е изд., перераб. и доп. Москва : Машиностроение,
1993. 640 с. ISBN 5-217-01304-0.
26. Вандевурд Д., Джосаттис Н. М., Д. Г. Шаблоны C++. Справочник раз­
работчика, 2-е изд. Санкт-Петербург : Альфа-книга, 2018. 848 с. ISBN
978-5-9500296-8-4.
27. Вейник А. И. Приближенный расчет процессов теплопроводности.
Москва : Госэнергоиздат, 1959. 184 с.
28. Влияние конфигурации и формы внешних ребер герметичных корпусов
технических средств на эффективность отведения тепла от процессора /
П. Г. Адамович [и др.] // Известия вузов России. Радиоэлектроника. 2023.
Т. 26, № 5. С. 63—75.
29. Галанин М. П., Родин А. С. Решение связанной задачи о термомеханиче­
ском контакте элементов твэла // Прикладная механика и техническая
физика. 2024. Т. 65, № 2. С. 99—109. DOI: 10.15372/PMTF202315387.
 95

30. Гусев А. А. Метод конечных элементов высокого порядка точности

решения краевых задач для эллиптического уравнения в частных про­
изводных. // Вестник РУДН. Серия МИФ. 2017. Т. 25, № 3. С. 217—233.
DOI: 10.22363/2312-9735-2017-25-3-217-233.
31. Деммель Д. Вычислительная линейная алгебра. Теория и приложения.
Москва : Мир, 2001. 430 с. ISBN 5-03-003402-1.
32. Донской А. А., Баритко Н. В. Кремнийорганические эластомерные теп­
лозащитные материалы низкой плотности // Каучук и Резина. 2003. № 2.
URL: [Link]
33. Евстафьев В. А. Конструирование космических аппаратов. Ч. 1: учебное
пособие. Санкт-Петербург : Балтийский государственный технический
университет, 2018. 99 с.
34. Зарубин В. С., Кувыркин Г. Н. Математические модели механики и элек­
тродинамики сплошной среды. Москва : Издательство МГТУ им. Н.Э.
Баумана, 2008. 512 с. ISBN 978-5-7038-3162-5.
35. Классман Е. Ю., Лутфуллин Р. Я. Влияние температуры нагрева заго­
товки перед теплой прокаткой на структуру и своства титанового сплава
BT22 // Фундаментальные проблемы современного материаловедения.
2024. Т. 21, № 2. С. 1811—1416. DOI: 10.25712/ASTU.1811-1416.2024.02.008.
36. Краснов М. М. Метапрограммирование шаблонов C++ в задачах мате­
матической физики. Москва : ИПМ им. М.В. Келдыша, 2017. DOI: 10.
20948/mono-2017-krasnov.
37. Кувыркин Г. Н. Математическая модель нелокальной термовязкоупругой
среды. Ч. 1. Определяющие уравнения // Вестник МГТУ им. Н.Э. Бау­
мана. Сер. Естественные науки. 2013. Т. 48, № 1. URL: [Link]
[Link]/articles/13/[Link].
38. Кувыркин Г. Н. Математическая модель нелокальной термовязкоупругой
среды. Ч. 2. Уравнение теплопроводности // Вестник МГТУ им. Н.Э. Ба­
 96

умана. Сер. Естественные науки. 2013. Т. 49, № 2. URL: [Link]

[Link]/articles/30/[Link].
39. Кувыркин Г. Н. Математическая модель нелокальной термовязкоупругой
среды. Ч. 3. Уравнения движения // Вестник МГТУ им. Н.Э. Баумана.
Сер. Естественные науки. 2013. Т. 50, № 3. URL: [Link]
ru/articles/119/[Link].
40. Кувыркин Г. Н., Савельева И. Ю. Численное решение интегро-диффе­
ренциального уравнения теплопроводности для нелокальной среды //
Математическое моделирование. 2013. Т. 25, № 5. С. 99—108.
41. Кувыркин Г. Н., Соколов А. А. Принцип Сен-Венана в задачах нело­
кальной теории упругости // Вестник МГТУ им. Н.Э. Баумана. Сер.
Естественные науки. 2023. Т. 109, № 4. С. 4—17. DOI: 10 . 18698 / 1812 -
3368-2023-4-4-17.
42. Кувыркин Г. Н., Соколов А. А. Решение задачи о напряженно-деформиро­
ванном состоянии пластины с эллиптическим вырезом при механических
и температурных нагружениях в нелокальной постановке // Прикладная
механика и техническая физика. 2024. № 4. С. 193—203. DOI: 10.15372/
PMTF202315385.
43. Лисовенко Д. С. Ауксетическая механика изотропных материалов, кри­
сталлов и анизотропных композитов : дис. ... д-ра физ.-мат. наук :
01.02.04. Москва: Иститут проблем механики им. А. Ю. Ишлинского Рос­
сийской Академии Наук, 2019. 392 с.
44. Лычев С. А. Законы сохранения недиссипативной микроморфной термо­
упругости // Вестник СамГУ. Естественнонаучная серия. 2007. Т. 54, № 4.
С. 225—262.
45. Морозов Н. Ф. Избранные двумерные задачи теории упругости. Ленин­
град : Издательство Ленинградского университета, 1978. 182 с.
 97

46. Печинкин А. В., Тескин О. И., Цветкова Г. М. Теория вероятностей :

учебник для втузов. Москва : Издательство МГТУ им. Н. Э. Баумана,
2006. 455 с. ISBN 5-7038-2485-0.
47. Писсанецки С. Технология разреженных матриц. Москва : Мир, 1988.
410 с. ISBN 5-03-000960-4.
48. Применение альтернативных серендиповых моделей при решении задач
о кручении призматических стержней. / И. А. Астионенко [и др.] // Вест­
ник ХНТУ. 2013. Т. 46, № 1. С. 356—361.
49. Савелов А. А. Плоские кривые: Систематика, свойства, примене­
ния. Справочное руководство. Москва : URSS, 2020. 294 с. ISBN
978-5-397-07388-2.
50. Савельева И. Ю. Вариационная формулировка математической моде­
ли процесса стационарной теплопроводности с учетом пространственной
нелокальности // Вестник МГТУ им.Н.Э.Баумана. Естественные науки.
2022. № 2. С. 68—86.
51. Савельева И. Ю. Влияние нелокальности среды на распределения тем­
пературы и напряжений в упругом теле при импульсном нагреве //
Известия РАН. Механика твердого тела. 2018. № 3. С. 45—52.
52. Савельева И. Ю. Двойственная вариационная модель стационарного
процесса теплопроводности, учитывающая пространственную нелокаль­
ность // Вестник МГТУ им.Н.Э.Баумана. Естественные науки. 2022. № 5.
С. 45—61.
53. Савельева И. Ю. Разработка и анализ математических моделей термоме­
ханики структурно-чувствительных материалов : дис. ... д-ра физ.-мат.
наук : 1.2.2. Москва: Московский государственный технический универси­
тет имени Н.Э. Баумана (национальный исследовательский университет),
2023. 375 с.
 98

54. Савельева И. Ю. Численное моделирование термоудара в упругом теле с

учетом эффектов нелокальности среды // Вестник МГТУ им. Н.Э. Бау­
мана. Естественные науки. 2020. № 3. С. 20—29.
55. Савин Г. Н. Распределение напряжений около отверстий. Киев : Наукова
Думка, 1968. 890 с.
56. Свидетельство о гос. регистрации программы для ЭВМ. NonLocFEM /
А. А. Соколов ; А. А. Соколов, И. Ю. Савельева. № 2021661966 ; заявл.
20.07.2021 ; опубл. 22.09.2022, РД040930 (Рос. Федерация).
57. Северюхин А. В., Северюхина О. Ю., Вахрушев А. В. Расчет ко­
эффициента теплопроводности нанокристаллов // Вестник Пермского
национального исследовательского политехнического университета. Ме­
ханика. 2022. № 1. С. 115—122. DOI: 10.15593/[Link]/2022.1.10.
58. Сен-Венан Б. Мемуар о кручении призм. Мемуар об изгибе призм. /
под ред. Г. Джанелидзе. Москва : Физматлит, 1961. 519 с.
59. Тамразян А. Г., Черник В. И. Жесткость поврежденной пожаром железо­
бетонной колонны при разгрузке после высокоинтенсивного горизонталь­
ного воздействия // Вестник МГСУ. 2023. Т. 18, № 9. С. 1369—1382. DOI:
10.22227/1997-0935.2023.9.1369-1382.
60. Трубицын В. Ю., Долгушева Е. Б. Особенности решеточной теплопровод­
ности наноструктурированных материалов на основе титана и алюминия.
Метод молекулярной динамики // Химическая физика и мезоскопия.
2019. Т. 21, № 4. С. 541—550. DOI: 10.15350/17270529.2019.4.57.
61. Численное моделирование задач термоупругости для конструкции с внут­
ренним источником / М. В. Васильева [и др.] // Математические заметки
СВФУ. 2018. Т. 24, № 3. С. 52—64. DOI: 10.25587/SVFU.2018.3.10889.

62. A novel stochastic photo-thermoelasticity model according to a diffusion in­

teraction processes of excited semiconductor medium / K. Lotfy [et al.] //
 99

The European Physical Journal Plus. 2022. Vol. 137. P. 721—738. DOI:
10.1140/epjp/s13360-022-03185-6.

63. Abdollahi R., Boroomand B. Benchmarks in nonlocal elasticity defined by

Eringen’s integral model // International Journal of Solids and Structures.
2013. Vol. 50, No. 18. P. 2758—2771. DOI: 10.1016/[Link].2013.04.027.

64. Abdollahi R., Boroomand B. Nonlocal elasticity defined by Eringen’s integral

model: Introduction of a boundary layer method // International Journal of
Solids and Structures. 2014. Vol. 51, No. 9. P. 1758—1780. DOI: 10.1016/j.
ijsolstr.2014.01.016.

65. Abdollahi R., Boroomand B. On using mesh-based and mesh-free methods

in problems defined by Eringen’s non-local integral model: issues and reme­
dies // An International Journal of Theoretical and Applied Mechanics. 2019.
Sept. Vol. 54. P. 1801—1822. DOI: 10.1007/s11012-019-01048-6.

66. Ahmadi G., Firoozbakhsh K. First strain gradient theory of thermoelastic­

ity // International Journal of Solids and Structures. 1975. Vol. 11, No. 3.
P. 339—345. DOI: 10.1016/0020-7683(75)90073-6.

67. Aifantis E. On the role of gradients in the localization of deformation and frac­
ture // International Journal of Engineering Science. 1992. Vol. 30, No. 10.
P. 1279—1299. DOI: 10.1016/0020-7225(92)90141-3.

68. Altan B. S., Aifantis E. C. On some aspects in the special theory of gradient
elasticity // Journal of Mechanical Behavior of Materials. 1997. Vol. 8, No. 3.
P. 31—282. DOI: 10.1515/JMBM.1997.8.3.231.

69. Altan S. B. Existence in nonlocal elasticity // Archive Mechanics. 1989.

Vol. 41. P. 25—36.

70. Altan S. B. Uniqueness in nonlocal thermoelasticity // Journal of Thermal

Stresses. 1991. Vol. 14. P. 121—128.
 100

71. Altan S. B. Uniqueness of initial-boundary value problems in nonlocal elastic­

ity // International Journal of Solids and Structures. 1989. Vol. 25, No. 11.
P. 1271—1278. DOI: 10.1016/0020-7683(89)90091-7.

72. Bathe K.-J. Finite Element Procedures. Second edition. 2014. 1065 p. URL:
[Link]
Edition_7th_Printing%20_1-[Link].

73. Bentley J. L. Multidimensional binary search trees used for associative search­
ing // Communications of the Association for Computing Machinery. New
York, USA, 1975. Vol. 18, No. 9. P. 509—517. DOI: 10.1145/361002.361007.

74. Benvenuti E., Tralli A. The fast Gauss transform for non-local integral FE
models // Communications in Numerical Methods in Engineering. 2006.
Vol. 22. P. 505—533. DOI: 10.1002/cnm.827.

75. Bogy D. B., Sternberg E. The effect of couple-stresses on the corner singularity
due to an asymmetric shear loading // International Journal of Solids and
Structures. 1968. Vol. 4, No. 2. P. 159—174.

76. Breakdown of Fourier’s Law in Nanotube Thermal Conductors / C. W. Chang

[et al.] // Physical review letters. 2008. Sept. Vol. 101. DOI: 10 . 1103 /
PhysRevLett.101.075903.

77. Cosserat E., Cosserat F. Theory of Deformable Bodies // A. Hermann et

Fils. 1909. P. 226.

78. Cuthill E., McKee J. Reducing the bandwidth of sparse symmetric matrices.
1969. DOI: 10.1145/800195.805928.

79. Determining the Elasticity of Materials Employing Quantum-mechanical Ap­

proaches: From the Electronic Ground State to the Limits of Materials
Stability / M. Friák [et al.] // steel research int. 2011. Vol. 82. P. 86—100.
DOI: 10.1002/srin.201000264.
 101

80. Duczek S. Higher order finite elements and the fictitious domain concept
for wave propagation analysis. Magdeburg, Germany : Otto von Guericke
University Library, 2014. P. 458. DOI: 10.25673/4151.

81. Edelen D. G. B., Green A. E., Laws N. Nonlocal continuum mechanics //

Archive for Rational Mechanics and Analysis. 1971. Vol. 43. P. 36—44.
DOI: 10.1007/BF00251544.

82. Edelen D. G. B., Laws N. On the thermodynamics of systems with nonlocal­

ity // Archive for Rational Mechanics and Analysis. 1971. Vol. 43. P. 24—35.
DOI: 10.1007/BF00251543.

83. Eremeyev V. A., Lazar M. Strong ellipticity within the Toupin–Mindlin first
strain gradient elasticity theory // Mechanics Research Communications.
2022. Vol. 124. DOI: 10.1016/[Link].2022.103944.

84. Eringen A. C. Linear theory of nonlocal elasticity and dispersion of plane

waves // International Journal of Engineering Science. 1972. Vol. 10, No. 5.
P. 425—435. DOI: 10.1016/0020-7225(72)90050-X.

85. Eringen A. C. Mechanics of micromorphic materials // Applied Mechanics /

ed. by H. Görtler. Berlin, Heidelberg : Springer Berlin Heidelberg, 1966.
P. 131—138. ISBN 978-3-662-29364-5. DOI: 10.1007/978-3-662-29364-5_12.

86. Eringen A. C. Microcontinuum field theories: foundations and solids.

NewYork : Springer-Verlag, 1999. 325 p. ISBN 978-0-387-98620-3. DOI:
10.1007/978-1-4612-0555-5.

87. Eringen A. C. Nonlocal continuum field teories. New York-Berlin-Heidelberg :

Springer-Verlag, 2002. P. 376. ISBN 978-0-387-95275-8. DOI: 10 . 1007 /
b97697.

88. Eringen A. C. Simple microfluids // International Journal of Engineering

Science. 1964. Vol. 2, No. 2. P. 205—217. DOI: 10.1016/0020- 7225(64)
90005-9.
 102

89. Eringen A. C., Edelen D. G. B. On nonlocal elasticity // International Journal

of Engineering Science. 1972. Vol. 10, No. 3. P. 233—248. DOI: 10.1016/
0020-7225(72)90039-0.

90. Eringen’s nonlocal and modified couple stress theories applied to vibrating ro­
tating nanobeams with temperature effects / A. Rahmani [et al.] // Mechanics
of Advanced Materials and Structures. 2021. Vol. 29, No. 26. P. 4813—4838.
DOI: 10.1080/15376494.2021.1939468.

91. Erosion behaviour of platinum aluminide bond coat on directionally solidified

CM247 and AM1 single crystal superalloys / S. L. Gokul [et al.] // Surface and
Coatings Technology. 2022. Vol. 429. DOI: 10.1016/[Link].2021.127941.

92. Fazilati J., Khalafi V., Shahverdi H. Three-dimensional aero-thermo-elasticity

analysis of functionally graded cylindrical shell panels // Proceedings of the
Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineer­
ing. 2019. Vol. 233, No. 5. DOI: 10.1177/0954410018763861.

93. Flatten A. Lokale und nicht-lokale Modellierung und Simulation thermo­

mechanischer Lokalisierung mit Schädigung für metallische Werkstoffe unter
Hochgescwindigkeitsbeanspruchungen. Berlin : der Bundersanstalt fur Mate­
rialforschung, 2008. P. 199.

94. Flaw Insensitive Fracture in Nanocrystalline Graphene / T. Zhang [et al.] //

Nano letters. 2012. Aug. Vol. 12. P. 4605—10. DOI: 10.1021/nl301908b.

95. Gao J. An asymmetric theory of nonlocal elasticity–Part 2. Continuum

field // International Journal of Solids and Structures. 1999. Vol. 36, No. 20.
P. 2959—2971. DOI: 10.1016/S0020-7683(97)00322-3.

96. Günther W. Zur Statik und Kinematik des Cosseratschen Kontinuums //

Abhandlungen der Braunschweigischen Wissenschaftlichen Gesellschaft Band.
1958. Vol. 10. P. 195—213. DOI: 10.24355/dbbs.084-201212211408-0.
 103

97. High Temperature Resistant Coatings for Strategic Aero space Applications /
Z. Alam [et al.] // Defence Science Journal. 2023. Vol. 73, No. 2. P. 171—181.
DOI: 10.14429/dsj.73.18638.

98. Hirschfelder J., Eyring H., Topley B. Reactions Involving Hydrogen Molecules
and Atoms // The Journal of Chemical Physics. 1936. Vol. 4, No. 3.
P. 170—177. DOI: 10.1063/1.1749815.

99. Hsu Y. C., Wang W. J. Couple-stress effects near an interior hole of an infi­
nite elastic plane subjected to a concentated force // Journal of the Franklin
Institute. 1973. Vol. 295, No. 5. P. 411—421.

100. Impact of time after fire on post-fire seismic behavior of RC columns /

U. Demir [et al.] // Structures. 2020. Vol. 26. P. 537—548. DOI: 10 .
1016/[Link].2020.04.049.

101. Kröner E. Elasticity theory of materials with long range cohesive forces // In­
ternational Journal of Solids and Structures. 1967. Vol. 3, No. 5. P. 731—742.
DOI: 10.1016/0020-7683(67)90049-2.

102. Kumar S., Haque A., Gao H. Notch insensitive fracture in nanoscale thin
films // Applied Physics Letters. 2009. June. Vol. 94. P. 253104—253104.
DOI: 10.1063/1.3157276.

103. Kuvyrkin G. N., Savelyeva I. Y., Sokolov A. A. 2D nonlocal elasticity: In­

vestigation of stress and strain fields in complex shape regions // Journal
of Applied Mathematics and Mechanics. 2023. Vol. 103, No. 3. DOI:
10.1002/zamm.202200308.

104. Kuvyrkin G. N., Savelyeva I. Y., Sokolov A. A. Features of the software

implementation of the numerical solution of stationary heat equation tak­
ing into account the effects of nonlocal finite element method // Journal of
Physics: Conference Series. 2020. Vol. 1479, No. 1. DOI: 10.1088/1742-
6596/1479/1/012034.
 104

105. Lazar M., Agiasofitou E. Toupin–Mindlin first strain-gradient elasticity for

cubic and isotropic materials at small scales // Proceedings in Applied Math­
ematics and Mechanics. 2023. Vol. 23. DOI: 10.1002/pamm.202300121.

106. Lei L., Pengfei L., Markus O. A molecular dynamics simulation study on
enhancement of thermal conductivity of bitumen by introduction of car­
bon nanotubes // Construction and Building Materials. 2022. Vol. 353.
P. 129—166. DOI: 10.1016/[Link].2022.129166.

107. Li X. The Coupling between Quantum Mechanics and Elasticity. the Faculty
of the Department of Mechanical Engineering University of Houston, 05/2016.
P. 110.

108. Madan G. K., Ronald L. H., Oswald B. Finite-element grid improvement

by minimization of stiffness matrix trace // Computers & Structures. 1989.
Vol. 31, No. 6. P. 891—896. DOI: 10.1016/0045-7949(89)90274-5.

109. Malvi B., Roy M. Elevated Temperature Erosion of Plasma Sprayed Thermal
Barrier Coating // Therm Spray Tech. 2021. Vol. 30. P. 1028—1037. DOI:
10.1007/s11666-021-01189-9.

110. Mathematical modeling of insulating coating of thermal conductivity includ­

ing body’s own radiation and non-local spatial effects / A. A. Sokolov [et al.] //
Journal of Physics: Conference Series. 2024. Vol. 2817, No. 1. P. 12—28.
DOI: 10.1088/1742-6596/2817/1/012028.

111. Maugis P. Nonlinear elastic behavior of iron-carbon alloys at the nanoscale //

Computational Materials Science. 2019. Vol. 152. P. 460—469. DOI: 10.
1016/[Link].2018.12.024.

112. Microstructure vs. Flaw: Mechanisms of Failure and Strength in Nanos­

tructures. / W. Gu [et al.] // Nano letters. 2013. Oct. Vol. 13. DOI:
10.1021/nl403453h.
 105

113. Mindlin R. D. Microstructure in Linear Elasticity // Archive for Rational Me­

chanics and Analysis. 1964. Vol. 16. P. 51—78. DOI: 10.1007/BF00248490.

114. Mindlin R. D. Second gradient of strain and surface-tension in linear elas­

ticity // International Journal of Solids and Structures. 1965. Vol. 1.
P. 417—438. DOI: 10.1016/0020-7683(65)90006-5.

115. Mindlin R. D. Stress function for a cosserat continuum // International Jour­

nal of Engineering Science. 1965. Vol. 1, No. 3. P. 265—271.

116. Mindlin R. D., Eshel N. N. On first strain-gradient theories in linear elas­

ticity // International Journal of Solids and Structures. 1968. Vol. 4.
P. 109—124. DOI: 10.1016/0020-7683(68)90036-X.

117. Mindlin R. D., Tierstin H. F. Effects of couple-stress in linear elasticity //

Experimental Mechanics. 1962. Vol. 11. P. 415—488.

118. Moosazadeh H., Mohammadi M. M. Time domain aero-thermo-elastic insta­

bility of two-dimensional non-linear curved panels with the effect of in-plane
load considered // SN Applied Sciences. 2020. Vol. 2. DOI: 10.1007/s42452-
020-03411-9.

119. Patrick M. K. Achieving finite element mesh quality via optimization of the
Jacobian matrix norm and associated quantities. Part I—a framework for
surface mesh optimization // International Journal for Numerical Methods
in Engineering. 2000. No. 48. P. 401—420. DOI: 10.1002/(SICI)1097-
0207(20000530)48:3<401::AID-NME880>[Link];2-D.

120. Pisano A. A., Sofi A., Fuschi P. Nonlocal integral elasticity: 2D finite ele­
ment based solutions // International Journal of Solids and Structures. 2009.
Vol. 46. P. 3836—3849. DOI: 10.1016/[Link].2009.07.009.
 106

121. Pisano A. A., Fuschi P., Polizzotto C. Integral and differential approaches
to Eringen’s nonlocal elasticity models accounting for boundary effects with
applications to beams in bending // Journal of Applied Mathematics and
Mechanics. 2021. Vol. 101, No. 8. DOI: 10.1002/zamm.202000152.

122. Polizzotto C. Nonlocal elasticity and related variational principles // Interna­

tional Journal of Solids and Structures. 2001. Vol. 38, No. 42. P. 7359—7380.
DOI: 10.1016/S0020-7683(01)00039-7.

123. Polizzotto C., Fuschi P., Pisano A. A. A nonhomogeneous nonlocal elasticity

model // European Journal of Mechanics - A/Solids. 2006. Vol. 25, No. 2.
P. 308—333. DOI: 10.1016/[Link].2005.09.007.

124. Rogula D. Introduction to Nonlocal Theory of Material Media //. Vienna :

Springer Vienna, 1982. P. 123—222. ISBN 978-3-7091-2890-9. DOI: 10 .
1007/978-3-7091-2890-9_3.

125. Ru C. Q., Aifantis E. C. A simple approach to solve boundary-value problems

in gradient elasticity // Acta Mechanica. 1993. Vol. 101. P. 59—68. DOI:
10.1007/BF01175597.

126. Structural Factors of Hardening of U8A Carbon Tool Steel under Cyclic Heat
Exposure / A. M. Guryev [et al.] // Technical Physics. 2023. Vol. 68, No. 8.
P. 171—176. DOI: 10.1134/S1063784223700020.

127. Sydnaoui I., Mohamed R., Ab Kadir M. A. Design Engineering The Effects
of Seasonal Thermal Stresses at Concrete Buildings in the Arabic Area //
Design Engineering (Toronto). 2022. Oct. Vol. 6. P. 721—738.

128. Temperature Stress Analysis of Super-Long Frame Structures Accounting for

Differences in the Linear Expansion Coefficients of Steel and Concrete / Y. Jia
[et al.] // Processes. 2021. Vol. 9, No. 9. DOI: 10.3390/pr9091519.
 107

129. Thermoelastic with photogenerated model of rotating microstretch semicon­

ductor medium under the influence of initial stress / A. Saeed [et al.] //
Results in Physics. 2021. Nov. Vol. 31. DOI: 10.1016/[Link].2021.104967.

130. Toupin R. A. Elastic materials with couple stresses. // Archive for Rational
Mechanics and Analysis. 1962. Vol. 11. P. 385—414. DOI: 10 . 1007 /
BF00253945.

131. Tuna M., Kirca M. Exact solution of Eringen’s nonlocal integral model for
bending of Euler–Bernoulli and Timoshenko beams // International Journal
of Engineering Science. 2016. Vol. 105. P. 80—92. DOI: 10.1016/[Link].
2016.05.001.

132. Türkmen İ., Yalamaç E. Effect of Alternative Boronizing Mixtures on Boride

Layer and Tribological Behaviour of Boronized SAE 1020 Steel // Metals and
Materials International. 2022. Vol. 28. P. 1—15. DOI: 10.1007/s12540-021-
00987-8.

133. Wang J., Altan S. B. Uniqueness in generalized nonlocal thermoelasticity //

Journal of Thermal Stresses. 1993. Vol. 16. P. 71—78.

134. Wen P., Huang X., Aliabadi M. Two Dimensional Nonlocal Elasticity Analysis
by Local Integral Equation Method // Computer Modeling in Engineering and
Sciences. 2013. Vol. 96, No. 3. P. 199—225. DOI: 10.3970/cmes.2013.096.
199.

135. Xianqiao W., James D. L. Micromorphic theory: a gateway to nano world //

International Journal of Smart and Nano Materials. 2010. Vol. 1, No. 2.
P. 115—135. DOI: 10.1080/19475411.2010.484207.

136. Yan Z., Cheng E. A Novel Monte Carlo Method to Calculate the Thermal
Conductivity in Nanoscale Thermoelectric Phononic Crystals Based on Uni­
versal Effective Medium Theory // Mathematics. 2023. Vol. 11, No. 5. DOI:
[Link]
 108

137. Yang N., Zhang G., Li B. Violation of Fourier’s law and anomalous heat
diffusion in silicon nanowires // Nano Today. 2010. Vol. 5, No. 2. P. 85—90.
DOI: 10.1016/[Link].2010.02.002.

138. Zienkiewicz O., Taylor R., Zhu J. Z. The Finite Element Method: Its Basis
and Fundamentals. Seventh edition. 2013. ISBN 978-1-85617-633-0. DOI:
10.1016/C2009-0-24909-9.
 109

Приложение

Ниже представлен Листинг 1 конфигурационного файла, на основе ко­

торого был произведён расчёт комбинированной задачи темплопроводности
и термоупругости из раздела 4.7. Конфигурационный файл выполнен в виде
структуры, описанной в формате JSONSchema, и содержит 6 основных полей:
task, save, mesh, thermal_boundaries, mechanical_boundaries и materials.
В поле task описаны основные характеристики запускаемой задачи: её
размерность, тип расчёта и зависимость расчёта от времени. В поле save содер­
жатся параметры сохранения результатов расчётов, указан путь сохранения
результатов, а также названия сохраняемых файлов и точность с которой
записывать результаты расчётов. В разделе mesh указан путь по которому нахо­
дится файл содержащий конечно-элементную сетку. Поля thermal_boundaries
и mechanical_boundaries содержат граничные условия для температурной и
механической задач соответственно. Здесь важно отметить, что граничные усло­
вия заданы на именованных границах, поэтому важно, чтобы сетка содержала
информацию об этих границах в виде групп элементов. В разделе materials ука­
заны физические параметры материала и модельные параметры отдельно для
уравнения теплопроводности и уравнения равновесия. Здесь также важно отме­
тить, что информация о материале задана на именованной группе элементов,
которые обозначают определённый материал. Таких групп может быть несколь­
ко, где для каждого материала могут быть заданы свои параметры.

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

ности и термоупругости
{
" task ": {
" dimension ": 2 ,
" problem ": " thermomechanical " ,
5 " time_dependency ": false
 110

},
" save ": {
" folder ": "/ path / to / save / folder " ,
" config ": " config " ,
10 " csv ": " solution_2d " ,
" vtk ": " solution_2d " ,
" precision ": 7
},
" mesh ": {
15 " path ": "/ path / to / mesh . su2 "
},
" thermal_boundaries ": {
" Left ": {
" kind ": " flux " ,
20 " flux ": -1
},
" Right ": {
" kind ": " flux " ,
" flux ": 1
25 }
},
" mechanical_boundaries ": {
" Left ": [
{ " pressure ": -1 } ,
30 null
],
" Right ": [
{ " pressure ": 1 } ,
null
35 ],
" Horizontal ": [
null ,
{ " displacement ": 0 }
]
40 },
 111

" materials ": {

" Material_Name ": {
" physical ": {
" conductivity ": 1.0 ,
45 " youngs_modulus ": 400.0 ,
" poissons_ratio ": 0.3 ,
" thermal_expansion ": 2.5 e -3
},
" thermal_model ": {
50 " local_weight ": 0.5 ,
" nonlocal_radius ": 0.2
},
" mechanical_model ": {
" local_weight ": 0.75 ,
55 " nonlocal_radius ": 0.1
}
}
}
}