WWW.DISSERS.RU

БЕСПЛАТНАЯ ЭЛЕКТРОННАЯ БИБЛИОТЕКА

   Добро пожаловать!

Pages:     || 2 | 3 |
-- [ Страница 1 ] --

Министерство образования и наук

и РФ Московский физико-технический институт (государственный университет)

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

Челноков Федор Борисович Численное моделирование деформационных динамических

процессов в средах со сложной структурой Специальность 05.13.18 Математическое моделирование, численные методы и комплексы программ ДИССЕРТАЦИЯ на соискание ученой степени кандидата физико-математических наук

Научный консультант: д. ф.-м. н., проф. Петров И. Б.

Москва – 2005 Оглавление Введение............................... 6 7 10 10 0.1. Способы дискретизации уравнений механики......... 0.2. Способы построения сетки в области интегрирования.... 0.2.1. Квадратная регулярная сетка............... 0.2.2. Регулярная сетка из четырехугольников, сопряженная с границами области интегрирования............ 0.2.3. Нерегулярная треугольная сетка............. 0.2.4. Изменение сетки при деформировании тел........ Глава 1. Численное решение уравнений упругости...... 1.1. Математическая модель..................... 1.2. Выбор системы координат.................... 1.3. Обобщение записи дифференциальных уравнений...... 1.4. Спектральное исследование системы.............. 1.4.1. Прямая задача....................... 1.4.2. Сопряженная задача.................... 1.4.3. Нормировка собственных векторов............ 1.4.4. Нулевые собственные значения.............. 1.4.5. Матрицы,, 1.................... 1.5. Покоординатное расщепление.................. 1.6. Разностные схемы......................... 1.7. Сеточно-характеристические схемы............... 1.8. Расчет на границе области интегрирования.......... 1.8.1. Заданная внешняя сила..................

12 14 14 18 19 20 22 23 24 26 29 29 30 31 35 39 42 Оглавление 1.8.2. Заданная скорость границы................ 1.8.3. Смешанные условия.................... 1.8.4. Условия поглощения и симметрии............ 1.8.5. Решение на границе при наличии правой части..... 1.9. Контакт между двумя телами.................. 1.9.1. Полное слипание...................... 1.9.2. Свободное скольжение................... 1.10. Интегрирование уравнений акустики.............. 1.11. Двумерные уравнения упругости................ 1.12. Эйлерова сетка и границы из маркеров............

47 48 49 52 54 54 55 56 57 58 61 63 63 65 68 74 77 79 81 85 89 91 92 Глава 2. Построение нерегулярной треугольной сетки.... 2.1. Представление триангуляции в программе........... 2.1.1. Наиболее компактный формат.............. 2.1.2. Расширенные структуры данных для ускорения поиска 2.2. Триангуляция невыпуклого многоугольника с полостями.. 2.3. Оптимальная триангуляция Делоне............... 2.4. Поддержание заданной плотности сетки............ 2.4.1. Сокращение граничных вершин.............. 2.5. Обоснование корректности алгоритма............. 2.6. Размеры внутренних треугольников сетки........... 2.7. Допустимые размеры длины граничного ребра........ 2.7.1. Минимальный угол границы тела............. 2.7.2. Обработка треугольников с двумя граничными ребрами 2.8. Трудоемкость поиска треугольника............... 2.9. Момент вырождения сетки при движении вершин с постоянной скоростью........................... 2.10. Примеры работы алгоритмов.................. Глава 3. Контакт между телами в динамических задачах 97. Оглавление 3.1. Поиск сегментов контактирующих границ........... 104 3.1.1. Структуры многомерного поиска............. 105 3.1.2. Триангуляция пространства между телами....... 106 3.2. Проверка сблизившихся узлов.................. 109 3.3. Несовпадение узлов в контактирующих телах......... 110 3.4. Примеры расчетов с контактом нескольких тел........ 113 Глава 4. Интерполяция в треугольнике............. 121 4.1. Реконструкция полинома заданного порядка......... 122 4.2. Кусочно-линейная интерполяция................ 127 4.3. Градиент интерполяционного полинома............ 130 4.4. Вычисление интеграла полинома в треугольнике....... 132 4.5. Монотонная квадратичная реконструкция........... 134 4.6. Борьба с ростом вариации при интерполяции......... 138 4.7. Сравнение численных методов, использующих регулярную и нерегулярную сетки........................ 140 4.7.1. Выполнение законов сохранения импульса и энергии.. 147 Глава 5. Численный метод для бесструктурных сеток.... 149 5.1. Уравнение переноса........................ 149 5.2. Гиперболическая система уравнений.............. 157 5.3. Сравнение одномерных схем на решении уравнения переноса Глава 6. Распространение упругих волн в массивных породах 167 6.1. Введение.............................. 167 6.2. Начальное состояние среды................... 171 6.3. Граничные условия........................ 173 6.3.1. Поверхности трещин.................... 174 6.4. Примененная неравномерная треугольная сетка........ 176 6.5. Исследование энергии в области интегрирования....... 177 6.6. Равномерность распределения полостей............ Оглавление 6.7. Оценка вариации плотности тела со случайным распределением круговых полостей..................... 181 6.7.1. Полости одного размера.................. 182 6.7.2. Случайное распределение размеров полостей...... 183 6.8. Детали численных экспериментов................ 185 6.9. Анализ результатов расчетов.................. 187 Глава 7. Распространение волн в перфорированных средах 7.1. Двумерная постановка задачи.................. 199 7.2. Трехмерная постановка задачи................. 201 Глава 8. Высокоскоростной удар по многослойной преграде 8.1. Постановка задачи........................ 211 8.2. Изменение скорости и положения шарика со временем.... 214 8.3. Подвижная расчетная сетка................... 215 8.3.1. Перестройка сетки как задача оптимизации....... 219 8.4. Учет разрушения материалов.................. 222 8.4.1. Результаты расчетов.................... 226 8.4.2. Увеличение рассчитываемого периода соударения за счет фрагментации........................ 228 Заключение.............................. 240 Список использованных источников.............. Введение В механике деформируемого твердого тела к настоящему времени разработано большое количество моделей [1–5], описывающих поведение сплошных сред, фазовые переходы в них, критерии разрушения и фрагментации тел под действием интенсивной нагрузки, а также континуальные модели развития разрушений. Часть этих моделей хорошо исследована и не ставится под сомнение, однако вычислить аналитическим образом вытекающие из них следствия можно лишь для бесконечно-малых воздействий и очень простых по форме тел, желательно обладающих различного рода симметриями. Исследование поведения реальных тел со сложной геометрией, подвергающихся значительным внешним воздействиям, приводящим к конечным деформациям, на основании этих моделей выполнить невозможно без привлечения компьютера. Другая часть моделей призвана описать наблюдаемые явления, такие как формирование сложного вида разрушений возле места сильного удара. Однако ответить на вопрос, действительно ли модель адекватно описывает данное явление невозможно, без проведения ряда компьютерных моделирований, называемых численными экспериментами. К задачам численного исследования стоит отнести и определение зачастую многочисленных параметров той или иной модели, которые практически невозможно измерить напрямую. Успешные попытки выполнения численного моделирования предпринимаются уже почти целый век, и за этот период достигнут существенный Введение прогресс в способах построения эффективных программ численного расчета. Однако часть трудностей, по-видимому, носит фундаментальный характер, и их преодоление актуально и сейчас. К ним можно отнести численное исследование многомерных динамических задач, включающее в качестве подзадачи необходимость описания мгновенного состояния тел (введение надлежащей сетки, методы интерполяции и т.д.), подвергающихся значительным деформациям и фрагментации на отдельные части. Не менее важным вопросом остается достижение высокой точности решения за приемлемое время на имеющемся компьютерном оборудовании, а также борьба с нефизичными эффектами, которые не следуют из сформулированной физической модели, а возникают в результате приближенного характера замены законов механики в интегральной или дифференциальной формах конечными соотношениями. Краткие обзоры-классификации используемых численных методов и видов сеток можно найти в работах [5, 6].

0.1.

Способы дискретизации уравнений механики К настоящему времени разработаны несколько принципиальных ме тодов численного решения дифференциальных уравнений в частных производных. Большинство этих методов являются сеточными. Существует множество толкований термина сетка. Здесь мы будем придерживаться следующего определения. Сеткой (mesh) в двумерном пространстве называют планарный граф, вершины которого соответствуют выбранным точкам (узлам сетки). Ребра графа составляют сеточные линии, а области, ограниченные ребрами, также носят названия ячеек сетки. Численный метод описывает состояние среды, ассоциируя данные с некоторыми объектами сетки (узлами, ребрами, ячейками). Для сеточных численных методов существует понятие шаблона (stencil) Введение — множества узлов сетки (или сеток для нескольких моментов времени для динамических уравнений), через значения в которых аппроксимируются уравнения в данной окрестности. Метод конечных разностей (nite-dierence method) основывается на замене производных в поставленных дифференциальных уравнениях на конечные, недифференциальные комбинации (в простейшем случае — разности) значений, определенных в точках сетки. Приближенное решение соn стоит из значений Uij (в двумерном случае), представляющих поточечное приближение к решению в узлах сетки (xi, yj ) во время tn :

n Uij u(xi, yj, tn ).

Метод конечных объемов (nite-volume method), с другой стороны, основывается на приближении средних значений u в ячейках сетки n Uij yj+1/2 yj1/ xi+1/ u(x, y, tn )dxdy.

xi1/ Для дифференциальных уравнений, которые возникают из законов сохранения и потоков сохраняющейся величины, это очень естественный подход. Дифференциальные уравнения часто возникают из более фундаментальной интегральной формы, констатирующей, что интеграл от u по конечному объему меняется со временем только вследствие потоков через границу этого объема. Если объем совпадает с ячейкой сетки, то это приводит к формуле обновления среднего по ячейке исходя из численных значений потока через границы ячейки. В этом заключается преимущество над конечно-разностной формулировкой, поскольку численное метод получается консервативным, точно выполняя законы сохранения. Программы, реализующие как метод конечных разностей, так и меn тод конечных объемов, могут хранить рассчитанные параметры Uij в виде прямоугольной таблицы. В первом случае в таблицу попадут значения в некоторых точках тела, во втором — усредненные значения по ячейкам.

Введение Рис. 1. Используемые значения в методах конечных разностей (слева) и конечных объемов (справа), стрелками показаны потоки между ячейками.

Но алгоритмы численных методов, записанные как действия над элементами таблицы, в обоих случаях будут очень похожи (рис. 1) и могут даже совпадать для линейных уравнений. Метод конечных элементов (nite element method), примененный, например, в [7–9], относится к вариационным методам [10]. Решение в нем ищется среди конечного множества функций, определенных согласованно с текущей сеткой. Метод конечных элементов наиболее часто используется в случае нерегулярности сетки из-за простоты его обобщения на этот случай. Методы конечных объемов и конечных элементов высоких порядков аппроксимации с упрощенной, аналитической реконструкцией для неструктурированных сеток получили названия методов спектральных объемов (spectral volume method) [11] и спектральных элементов (spectral element method) [12] соответственно. Некоторые соображения, лежащие в их основе, были заимствованы в данной работе для построения сеточно-характеристического метода второго порядка на неструктурированных сетках. Методы частиц [13] и гладких частиц (smoothed particle hydrodynamics) [14–16], или метод свободных точек [17] применяются при моделировании значительных деформаций тел и разрушений в них, приводящих к отколу множества обломков. Сеточно-характеристический метод (grid-characteristic method) [18– 34], используемый в данной работе, ближе всего к конечно-разностным методам в том смысле, что в точках сетки хранятся приближенные зна Введение чения численного решения в этих точках, и в качестве основы используется дифференциальная форма записи уравнений механики. Однако вместо производных, входящих в исходные уравнения, аппроксимации подвергаются производные вдоль характеристических направлений, что обеспечивает большую точность. Также отличительной особенностью сеточнохарактеристического метода является удобство вычисления на границе области интегрирования с привлечением ровно того же числа граничных условий, которые необходимы для корректного замыкания системы дифференциальных уравнений.

0.2.

Способы построения сетки в области интегрирования Методы введения сеток для расчета деформируемых тел классифи цируются по трем основным признакам: • наличию связи между положением сеточных линий и границ исследуемых тел, • способностью меняться по мере деформирования моделируемых тел и • наличию структуры (регулярности) в сетке. Рассмотрим далее возможные способы введения сетки на примере плоского тела простейшей треугольной формы. 0.2.1. Квадратная регулярная сетка Рис. 2. Квадратная регулярная сетка в окрестности исследуемого объекта.

Введение Квадратная регулярная сетка (иногда употребляется термин декартова решетка (cartesian grid) — наиболее простой вид сетки, которая строится независимо от положения границ моделируемых тел (рис. 2). Достоинством этого подхода является исключительная простота вычисления положения точек сетки, в результате чего можно даже не хранить их координаты в программе, а вычислять каждый раз налету. Кроме того, одинаковая квадратная форма всех ячеек и равное расстояние между точками сетки значительно упрощает запись численных методов.

Рис. 3. Декартова решетка вблизи границы области интегрирования. Слева отмечен несимметричный шаблон конечно-разностных схем, справа — ячейка неправильной формы конечнообъемных методов.

К недостаткам решетки следует отнести сложность построения численных методов вблизи границ объектов. Здесь обычно избирается один из двух способов: использование точек или ячеек, не лежащих целиком внутри тела, либо введение дополнительных точек сетки в местах пересечения линий решетки с границами (рис. 3). В первом варианте приходится пользоваться каким-либо методом экстраполяции для определения значений в точках, лежащих за пределами рассчитываемого тела [6], что нелегко сопрячь с поставленными граничными условиями. Во втором варианте в конечно-разностном подходе расстояния до точек сетки вблизи границы могут сокращаться до бесконечно малых величин, а в конечно-объемном подходе также образовываться ячейки сложной формы. Поскольку размеры минимальной ячейки в задачах динамики диктуют ограничение на допустимый шаг интегрирования, то эта особенность становится особенно неприятной. В любом случае построение алгоритмов для расчета возле границы не является столь же простым делом, как для внутренних точек, Введение а также приводит к ухудшению качества решения в сравнении с другими подходами — утрачивается основное достоинство квадратных регулярных сеток. Поэтому такой метод применяют в основном, когда границы области интегрирования совпадают с осями выбранной декартовой системы координат. 0.2.2. Регулярная сетка из четырехугольников, сопряженная с границами области интегрирования Под регулярной сеткой (grid) обычно понимается такая сетка, через каждую точку которой проходит ровно две сеточные линии (у внутренних узлов — четыре примыкающих ребра), и существует взаимнооднозначное непрерывное отображение точек регулярной решетки на точки декартовой решетки.

Рис. 4. Возможные способы построения регулярной сетки в треугольной области.

Существуют алгоритмы [35, 36], позволяющие так разместить точки регулярной сетки (body-tted grid), чтобы они не выходили за пределы тела, а образованные ячейки между линиями сетки были максимально близки по размерам друг к другу и похожи по форме на квадрат. При этом сетка и шаблоны численных методов не отличаются принципиально внутри и на границе области интегрирования. На рис. 4 представлены возможные виды регулярной сетки для области треугольной формы. Предлагаемые алгоритмы обычно решают задачу оптимизации итерационными методами, последовательно приближаясь к оптимальному значению. Как следствие, полученная сетка существенно зависят от начального приближения, а при плохом выборе не всегда обеспечивается сходимость алгоритма Введение и отыскание глобального экстремума. Сложность построения регулярной сетки в телах сложной геометрии заставляет некоторых исследователей [6] возвращаться к декартовой решетке с ее упомянутыми недостатками. Даже для тела такой простой треугольной формы не могут одновременно гарантироваться следующие важные для качества решения параметры: • отсутствие острых углов между пересекающимися сеточными линиями, • малая кривизна сеточных линий, • примерно равные площади ячеек, расстояния между узлами вдоль сеточных линий и плотность узлов сетки во всех ее частях.

Рис. 5. Сетка, части которой являются регулярными.

На рис. 5 представлена еще одна сетка, через каждую точку которой проходит ровно две сеточные линии, однако взаимнооднозначное отображение на решетку отсутствует, и, следовательно, сетка не может считаться регулярной. Такого рода сетки возникают при разбиении исходной области на ряд четырехугольных частей, в каждой из которых сетки строятся независимо, после чего они «сшиваются». Если конечно-разностному численному методу для расчета узла важно лишь наличие четырех расположенных симметрично соседних узлов, то такого рода сетка может быть использована вместо регулярной. Однако, например, для методов, в которых уравнения «расщепляются» вдоль сеточных линий, эта сетка не подойдет.

Введение 0.2.3.

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

Рис. 6. Нерегулярная сетка в треугольной области.

Чаще всего используются нерегулярные сетки (рис. 6) с простейшей формой ячейки — треугольником, хотя иногда применяют четырехугольные и смешанные (треугольно-четырехугольные) нерегулярные сетки. Ячейки более сложной формы используются намного реже. 0.2.4. Изменение сетки при деформировании тел Для решения эволюционных задач выбираются несколько промежуточных моментов времени, называемых временными слоями, между исходным состоянием и полным временем моделирования. Решение последовательно отыскивается на каждом слое, причем расчетная сетка может также меняться. В совокупности сетки из всех слоев образуют единую сетку в пространстве координат-времени. Численные методы также классифицируются по виду движения узлов между отдельными сетками при переходе Введение к следующему слою. В механике сплошной среды принято записывать динамические уравнения в одном из двух видов переменных [5]: эйлеровых — фиксированных в пространстве и лагранжевых — движущихся вместе с точками тела. Аналогично выделяют и два основных типа сеток, но кроме них используют и другие [19]. Рассмотрим далее каждый из них отдельно. • Эйлерова сетка строится однократно и в дальнейшем не изменяется. Это может быть наилучшим решением для тел с фиксированными границами. Для решения задач с конечными деформациями границ выбирают наиболее простую форму эйлеровой сетки — декартову решетку, потому что, независимо от формы ячеек, неподвижные узлы и ребра сетки не будут совпадать в каждый момент времени с движущимися границами, и неизбежны сложности, описанные для декартовых решеток. • Точки лагранжевой сетки смещаются вместе с точками тела, поэтому границы тела всегда совпадают с сеточными линиями. Однако при наличии сдвиговых деформаций в теле, ячейки лагранжевой сетки могут постепенно вырождаться и пересекать друг друга. Численные методы при вырождении сетки обнаруживают неустойчивость и счет приходится прекращать. Поэтому неизбежен подход, называемый лагранжева сетка с перестройкой, в котором, как только детектируется приближение сетки к вырожденному состоянию, производится построение новой сетки, а значения в новых узлах интерполируются из прежних значений. Процедура интерполяции сама по себе приводит к потере точности решения, поэтому желательно перестраивать сетку как можно реже. В противном случае имеет смысл обратиться к подвижной сетке. Задачи о соударении в лагранжевых переменных рассматривались в работах [37, 38].

Введение • Подвижная сетка является обобщением лагранжевой, ее точки движутся от слоя к слою с ненулевыми скоростями и смещаются как относительно неподвижной системы координат, так и относительно точек тела. Их движение может быть подобрано таким образом, чтобы исключить вырождение сетки со временем. Однако использование подвижной сетки требует вносить изменения в решаемые уравнения для учета конвективных членов. • По-координатная комбинация разных типов движения сетки рассматривается в статье [19]. Например, если правая и левая границы области интегрирования неподвижны, то x-координата узлов может быть эйлеровой и не храниться в памяти программы, а y-координата — лагранжевой. • Совместная эйлерово-лагранжева сетка представлена в работе [39]. Подобная сетка имеет неподвижные точки внутри тела и движущиеся вместе с поверхностью тела точки на ее границе. Основная сложность работы с такими сетками — обеспечить сопряжение двух типов узлов в приграничной полосе. В другой работе [40] исследование задач о соударении в эйлеровых переменных приводило к необходимости использовать метода маркеров. • Иной вариант совмещения эйлерового и лагранжевого подхода предлагается в методе «частиц в ячейке» Ф. Харлоу [13,17]. Метод использует две сетки частиц: лагранжеву сетку для описания частиц, представляющих элементы жидкости, движущихся по неподвижной эйлеровой сетке. Другой пример совмещения неподвижной сетки с подвижными частицами дает работа [41]. • Неструктурированная лагранжева сетка с перестройкой (изменением «соседства») используется в методике «Медуза», также изложенной в [17]. Область интегрирования разбивается на построенные вокруг Введение заданного набора точек многоугольные ячейки («глобулы»), во многом схожие с ячейками диаграммы Вороного [42]. Заметную трудность в этой методике представляет расчет граничных ячеек, поскольку их стороны не совпадают с контуром области интегрирования. • В данной работе предлагается численный метод, основанный на нерегулярной лагранжевой сетке с локальной перестройкой. Для достижения порядков аппроксимации выше первого в этом методе рассчитываются узлы, находящиеся не только в вершинах треугольников сетки, но также и дополнительные узлы, например, лежащие в центре ребер. Поэтому первые узлы являются лагранжевыми, а вторые — подвижными.

Глава 1 Численное решение уравнений теории линейной упругости Численное решение уравнений линейной теории упругости представляет неподдельный интерес для механики деформируемого твердого тела, несмотря на то, что реальные материалы удовлетворяют этим уравнениям лишь при незначительных быстрых нагрузках. Дело в том, что численное исследование сложных реологических моделей (упругопластической, упруговязкопластической и других сред) применением метода расщепления по физическим процессам может быть сведено к некоторой корректировке упругого решения [43]. В настоящей главе рассматриваются сеточно-характеристические схемы решения уравнений теории упругости в случае произвольного выбора базиса в пространстве, который зависит от локальной формы разностной сетки. Представленная явная форма записи схем имеет вычислительные преимущества перед стандартной формулировкой [19]. Особое внимание уделено решению на поверхности области интегрирования при различных типах граничных условий, а также случаю контакта двух упругих сред.

Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ 1.1.

Математическая модель Сформулируем здесь уравнения, которым согласно [5] подчиняется состояние бесконечно малого объема сплошной линейно-упругой среды. Вопервых, это локальное уравнение движения: v = · T. Здесь — плотность материала, v — скорость движения среды в данной точке, — градиент по пространственным координатам, T — тензор напряжений Коши. В силу закона парности касательных напряжений [5] T является симметричным тензором. Компоненты свертки градиента с этим тензором задаются равенством ( · T )k = (T · )k = ( · T)k = i T T i Tki, где i — скалярный дифференциальный оператор. Здесь и далее предполагается, что индексы (i, j, k, l,... ) перебирают числа 1,2,3. В записи (T · ) необходимо помнить, что дифференциальный оператор действует на T, а не на последующие множители. Введем симметричный тензор малых деформаций 1 e = ( u + u ), 2 где u — поле перемещений (x = X + u), — оператор тензорного произведения: (a b)ij = ai bj. В произведении (u ) дифференцированию также подвергается u, а не последующие множители: ( u + u )ij = i uj + j ui. (Ниже по тексту поле перемещений больше не рассматривается, а символ u будет использоваться в другом смысле.) Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Линейная упругость материала подразумевает следующую связь напряжений с деформациями (закон Гука): T = (e : I)I + 2µe, (1.1) где, µ — параметры Ляме, определяющие свойства упругого материала (существуют взаимооднозначные формулы связи параметров Ляме с одной стороны и модуля Юнга и коэффициента Пуассона с другой), I — единичный тензор, «:» — двойное скалярное произведение1 (свертка): a:b= ij aij bij.

Чтобы получить замкнутую систему уравнений относительно скоростей и напряжений остается только продифференцировать по времени связь напряжений с деформациями. v = · T, T = ( · v)I + µ( v + v ). (1.2) 1.2.

Выбор системы координат Особенностью численного интегрирования дифференциальных урав нений можно назвать тот факт, что в программной реализации бывает удобно хранить состояние среды (в данном случае это компоненты скорости и напряжения) в некоторой фиксированной декартовой системе координат (x1, x2, x3 ). В то же время производные приближенно вычисляются вдоль направлений (1, 2, 3 ), выбираемых в зависимости от локальной конфигурации узлов сетки. При расчете на регулярной кубической сетке такая система координат возникает, если ввести (не обязательно единичные) базисные векторы, исходя из координат шести ближайших к данному узлов r,i, r+,i : i = r+,i r,i. двойная свертка с единичным тензором — это сумма диагональных элементов тензора или его след.

Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ В этом же случае непосредственно вычислимы производные функции u, значения которой известны в узлах сетки, лишь вдоль направлений i : u+,i u,i u. i 2 Приближенные равенства справедливы с точностью O( i 2 ) в центре отрезка, соединяющего r,i и r+,i, и с точностью O( i ) в прочих точках. Чтобы найти производную u вдоль других направлений, рассмотрим равенства (x1, x2, x3 ) = 1 2 3, (1, 2, 3 ) (1, 2, 3 ) = (x1, x2, x3 ) qi · j = ij. Здесь было обозначено j. x То есть (q1, q2, q3 ) является биортогональным базисом по отношению к (1, 2, 3 ): qj (1, 2, 3 ) (1, 2, 3 ) (1, 2, 3 ) здесь в числителе имеет место векторное произведение, а в знаменателе — смешанное. Для заметного упрощения последующих формул свяжем c каждым из векторов qj ортонормированный базис (nj,0, nj,1, nj,2 ), в котором nj,0 сонаправлен с qj : nj,0 = а единичные вектора nj,1, nj,2 восстановить из контекста. Уравнения (1.2) записаны в форме, независимой от выбора той или иной системы координат. Необходимо лишь специально указать, какие операторы подставлять вместо компонент градиента. В декартовой системе qj, lj = qj, qj выбираются произвольно. Далее мы будем q1 = [2, 3 ], q2 = [3, 1 ], q3 = [1, 2 ], (x1, x2, x3 ) (1, 2, 3 ) T q1 T = q2, T q опускать первый индекс j, а в nj,0 — второй индекс 0, если потерю можно Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ координат = {1, 2, 3 } =, x i =. xi. x В ином случае применим правила замены переменных: i = j j, xi j = {1, 2, 3 } = (1.3) 1.3.

Обобщение записи дифференциальных уравнений В вычислительной математике для анализа системы уравнений и по строения разностных схем для нее принято приводить ее запись к канонической форме u+ j Aj u = 0, j (1.4) где в вектор u собираются все переменные системы, в нашем случае (1.2) это u = {v, T} = {v1, v2, v3, t11, t12, t13, t22, t23, t33 }. Чтобы получить запись матриц Aj, представим оператор градиента (1.3) в виде суммы трех операторов, каждый из которых содержит производные лишь по одному направлению j : = j T T j, j = j = qj = lj nj. x j j j j Подстановка = j lj nj (1.5) в (1.2) приводит к следующей записи действия матриц на произвольный вектор: 1 v (nj · T). Aj = lj (nj · v)I + µ(nj v + v nj ) T (1.6) В запись разностных схем входит как сама матрица Aj, так и различные функции от нее. Чтобы иметь возможность в явном виде вычислять Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ разностные схемы, исключая процедуры, например, обращения матриц и приближенного вычисления ее собственных величин, проводится следующее исследование свойств Aj. Желание получить явную и максимально упрощенную запись разностных схем мотивируется стремлением эффективного использования машинных ресурсов и оптимизацией программы с точки зрения размеров и скорости счета, а также повышением ее точности за счет уменьшения ошибок округления и даже надежности за счет исключения методов, которые могут завершиться остановкой программы, например, обращение вырожденной матрицы.

1.4.

Спектральное исследование системы Целью данного раздела является нахождение всех собственных знаT чений c и собственных векторов u = {v, T} матриц Aj : Au = cu (1.7) (до конца раздела нижний индекс j опускается). Собственные значения должны иметь размерность скорости, чтобы приведенное выражение было согласованным с точки зрения размерности. Если предположить, что у матриц A есть полный набор собственных значений {ci } и собственных векторов {ui } (что соответствует действительности, и будет показано ниже), то можно ввести следующие обозначения: = diag{ci } — диагональная матрица, составленная из собственных значений A, 1 = {ui } — матрица, столбцами которой являются собственные векторы, следующие в том же порядке, что и собственные значения в. Поскольку полный набор собственных векторов образует линейно независимый базис, то имеется Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ полное право обозначить невырожденную матрицу, образованную ими, как обратную от некоторой другой невырожденной матрицы. Условие того, что каждый столбец 1 является собственным вектором A можно записать как A1 = 1. После домножения справа на : A = 1, и домножения слева на : A =. Из последнего равенства видно, что каждая строка является собственной строкой A или, что то же самое, после транспонирования — собственным вектором A. Найти эти собственные векторы можно из решения сопряженной задачи по отношению к (1.7): A u = cu.

T T (1.8) Из приведенных выкладок следует, что сопряженная задача будет иметь тот же набор собственных значений, а собственные вектора можно подобрать согласованно: 1 = I. Отсюда: размерности элементов — обратные скорости и напряжения, которые после домножения на 1 дают безразмерные единицы в I. 1.4.1. Прямая задача Возвращаясь от обобщенной записи (1.7) через подстановку (1.6) к переменным v и T, получаем систему определяющих уравнений n · T = l1 cv, (n · v)I + µ(n v + v n) = l cT.

(1.9) Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ После подстановки одного уравнения в другое (n · v)n + µn · (n v + v n) = l2 c2 v и учета тождества n · (n v + v n) = (n · v)n + v приходим к ( + µ)(n · v)n = (l2 c2 µ)v. Последнее уравнение имеет три качественно различные группы решений. 1. Векторы v и n коллинеарны, тогда c2 = l2 +2µ. После введения обозначения (продольная скорость звука в упругом материале) c1 = + 2µ T = (I + 2µn n), c1 (1.10) получаем пару собственных значений и векторов c = ±lc1, v = n, где — произвольный ненулевой множитель с размерностью скорости. 2. Векторы v и n ортогональны и v = 0. Левая часть уравнения обращается в нуль. Тогда оно справедливо только лишь, когда правая часть обращается в нуль также: c2 = µ l2. После введения другого обозначения (поперечная скорость звука) c2 = µ (1.11) можно описать это семейство собственных значений и векторов c = ±lc2, n · v = 0, µ T = (n v + v n). c Уравнение n · v = 0 имеет два линейно независимых решения, следовательно, всего в данном семействе четыре линейно независимых Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ собственных вектора {v, T}, но лишь два собственных значения, отличающихся знаком. 3. Собственное значение равно нулю и компоненты собственного вектора, сопоставленные скорости равны нулю. Возвращаясь к (1.9) получаем уравнение на оставшиеся компоненты собственного вектора: c = 0, v = 0, n · T = 0. (1.12) Равенство n · T = 0 содержит три скалярных произведения векторов, учитывая, что в T шесть независимых компонент, получается три линейно независимых собственных вектора, отвечающих нулевому собственному значению. В пункте (1) найдено 2 собственных вектора, в пункте (2) — 4 собственных вектора и еще 3 — в пункте (3). Итого, найдены все 9 линейно независимых пар: собственное значение – собственный вектор. 1.4.2. Сопряженная задача Операция транспонирования матрицы Aj (1.6) определяется из равенства u1 · Aj u = u · Aj u1. Раскроем правую часть и перегруппируем множители в ней: v1 v 1 · Aj = lv · (nj · T1 ) + lT : [(nj · v1 )I + µ(nj v1 + v1 nj )] = T1 T 1 = lT1 : (v nj + nj v) + lv1 · [(T : I)nj + 2µnj · T]. 2 Следовательно v (T : I)nj + 2µnj · T T. (1.13) Aj = lj 1 T 2 (v nj + nj v) Как уже обсуждалось после (1.8), величины v и T здесь имеют размерности обратных скоростей и напряжений соответственно.

T Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Рассуждая аналогичным образом, приходим к трем группам собственных значений и собственных векторов: c = ±lc1, c = ±lc2, c = 0, v = n, T= n n;

c n · v = 0, v = 0, T= 1 (n v + v n);

2c2 (1.14) (T : I)n + 2µn · T = 0.

Важные обозначения и соотношения Ненадолго прервемся и посвятим данный подраздел тому, что введем дополнительные обозначения, а также сформулируем наиболее примечательные тождествам между введенными векторами и матрицами, на которые будем ссылаться многократно ниже по тексту. Определим симметричные матрицы Nij = Nji как 1 Nij = (ni nj + nj ni ). (i, j = 0, 1, 2) 2 Матрицы Nij ортогональны между собой относительно операции двойного скалярного произведения: i j, k l : Nij : Nkl = ik jl 1 + ij. Последние соотношения легко проверить, если обратиться, в частности, к тождествам a a : I = a 2, ab:ab= a b 2, a b : b a = (a · b)2, a a : a b = a a : b a = a 2 a · b, справедливым для любых векторов a и b. Нетрудно показать также, каково будет действие Nij на вектор nk : i, j, k : 1 Nij · nk = ((nj · nk )ni + (ni · nk )nj ). Следовательно N11, N12, N22, а также всевозможные их линейные комбинации при действии на n дают нуль, а значит, являются собственными векторами прямой задачи.

Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Другими важным свойствам матриц Nij является результат двойной свертки с произвольной симметричной матрицей T i, j : Nij : T = ni · T · nj, (1.15) что в случае единичной матрицы дает i, j : Nij : I = ni · nj = ij. (1.16) N00 + N11 + N22 = I.

Разложение произвольного вектора v по ортонормированному базису: v = (n · v)n + (n1 · v)n1 + (n2 · v)n2, приводит к 1 (n v + v n) = (n · v)N00 + (n1 · v)N01 + (n2 · v)N02. 2 (1.18) (1.17) Система N00, N01, N02, N12, N11 N22, N11 + N22 также ортого нальна, но не полностью нормирована ( N = N : N), следовательно, разложение произвольной симметричной матрицы T имеет вид T =(N00 : T)N00 + 2(N01 : T)N01 + 2(N02 : T)N02 + 2(N12 : T)N12 + 1 1 + ((N11 N22 ) : T)(N11 N22 ) + ((N11 + N22 ) : T)(N11 + N22 ). 2 2 (1.19) Применение (1.17) к вектору T · n дает T · n = (n · T · n)n + (n1 · T · n)n1 + (n2 · T · n)n2 = = (N00 : T)n + (N01 : T)n1 + (N02 : T)n2. Применение (1.18) к вектору T · n дает 1 ((T · n) n + n (T · n)) = (T : N00 )N00 + (T : N01 )N01 + (T : N02 )N02. 2 (1.21) (1.20) Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ 1.4.3.

Нормировка собственных векторов Скалярное произведение собственных векторов прямой задачи {vd, Td }, составляющих столбцы 1, и сопряженной задачи {vc, Tc }, составляющих строки, должно равняться единице для соответствующих пар векторов (и нулю для прочих пар): vd · vc + Td : Tc = 1. Подставим в это равенство найденные выше значения собственных векторов, чтобы согласовать их длину. Для первой группы собственных векторов имеем c = ±lc1, vc = n, Tc = 1 N00, c l2 1+l2.

vd = n, Td = (I + 2µN00 ), c где длина собственного вектора сопряженной задачи была произвольным образом фиксирована, а = Для второй группы — c = ±lc2, vc = nk, Tc = 1 N0k, c2 vd = nk, Td = 2 µN0k, c (k = 1, 2) при том же значении. 1.4.4. Нулевые собственные значения Вектора {0, Tc,i }, где Tc,1 = N12, Tc,2 = N11 N22, Tc,3 = N11 + N22 2 N00, + 2µ удовлетворяют условиям на собственный вектор сопряженной задачи (1.14). Система биортогональных к ним векторов, подходящих прямой задаче (1.12), такова: Td,1 = 2N12, 1 Td,2 = (N11 N22 ), 2 1 Td,3 = (N11 + N22 ). Коэффициенты подбирались так, чтобы i, j : Tc,i : Td,j = ij.

Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ 1.4.5.

Матрицы,, Можно подвести итог предшествующих подразделов, представив явный вид матриц, на которые раскладывается A: = l diag{c1, c1, c2, c2, c2, c2, 0, 0, 0}, 1 w1 n · v c1 N00 : T w 1 n · v + c1 N00 : T 2 w 1 n1 · v c2 N01 : T 3 w 1 n1 · v + c2 N01 : T 4 v 1 w = = n2 · v c2 N02 : T 5 T 1 w6 n2 · v + c2 N02 : T w7 N12 : T w8 (N11 N22 ) : T w9 (N11 + N 2 +2µ N00 ). (1.22) (1.23) :T К сожалению, исследование собственных векторов системы (1.7) позволяет выписать лишь действие матрицы 1 на вектор, находящийся слева: w 1, но не справа: 1 w, для чего потребуются дополнительные усилия. Чтобы определить u = 1 w, будем решать систему уравнений u = w. После сложения трех первых пар строк (1.23) n·v = w 1 + w2, 2 n1 · v = w 3 + w4, 2 n2 · v = w 5 + w6. T Подставляя это в (1.17), имеем v= w 5 + w6 w3 + w 4 w 1 + w2 n+ n1 + n2. 2 2 2 w4 w 3, 2 w 6 w5. Теперь попарно вычтем те же строки системы уравнений: N00 : T = c1 w 2 w1, 2 N01 : T = c2 N02 : T = c Последние строки дают N12 : T = w7, (N11 N22 ) : T = w8, (N11 +N22 ) : T = w9 + c1 (w2 w1 ). + 2µ Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Также подставим это в (1.19) и сразу заменим сумму N11 + N22 по формуле (1.16):

(w1 + w2 )n + (w3 + w4 )n1 + (w5 + w6 )n2 (w2 w1 )[(c1 c3 )N00 + c3 I]+ v 1. = 1 w = 2 +2c2 (w4 w3 )N01 + 2c2 (w6 w5 )N02 + T +4w7 N12 + w8 (N11 N22 ) + w9 (I N00 ) (1.24) В дальнейшем нам потребуются оперировать отдельными столбцами 1. Чтобы получить каждый из них, достаточно помножить матрицу на соответствующий единичный вектор:,i 1 ei, или, что то же самое, положить в записи (1.24) wj = ij. Постолбцовая запись 1 : n1 n n n1 1 1 = 2 [(c1 c3 )N00 + c3 I] [(c1 c3 )N00 + c3 I] 2c2 N01 2c2 N01 n2 n2 0 0 0. 2c2 N02 2c2 N02 4N12 N11 N22 I N00 (1.25) 1.5.

Покоординатное расщепление Запись уравнений в форме (1.4) удобна тем, что позволяется постро ить разностную схему интегрирования в трехмерном пространстве, если исследованы одномерные схемы для системы уравнений u = 0. u + Aj j следовательными временными слоями: un+1 = fj (Aj )un. (1.26) Обозначим действия одномерной схемы оператором перехода f между по Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ В дальнейшим будем отождествлять схему и соответствующий ей оператор f, например, говоря: «схема f ». Для всех одномерных явных схем, приведенных далее в данной работе, справедливо ограничение на максимальный шаг интегрирования по времени: j < jmax = 1 cmax j = 1, c1 lj где cmax — максимальное по модулю собственное значение Aj, c1 — проj дольная скорость звука (1.10). Также в этом условии предполагается, что шаг между узлами сетки в направлении j равен единице. Подход, состоящий в комбинировании операторов перехода для одномерных систем (1.26), чтобы построить разностные схемы для трехмерных систем (1.4), носит название расщепление по пространственным координатам. Приведем далее (без доказательства) несколько принципиальных способов построения схем расщеплением, каждый из которых использовался в последующих главах. Простейший способ Трехмерный оператор перехода (un+1 = F (A1, A2, A3 )un ) F (A1, A2, A3 ) = j j fj 1 Aj j (1.27) при условиях j = 1, j j > обеспечивает первый порядок аппроксимации, если в качестве базисных (fj ) были взяты схемы как минимум первого порядка. Набор {j } подбирается из условия максимизации допустимого шага интегрирования cmax cmax cmax 1 =2 =3, 1 2 3 которое приводит к cmax i. i = max c1 + cmax + cmax 2 Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ и < max = j cmax j = c j lj.

(1.28) Видно, что максимально допустимый шаг интегрирования, тем не менее, снижается по сравнению с одномерным случаем. Если операторы fj (Aj ) являются линейными функциями своего аргумента, что, как правило, имеет место для схем первого порядка точности, то (1.27) заметно упрощается: F (A1, A2, A3 ) = j fj (Aj ).

Хотя величины {j } уже не фигурируют в формуле, но более строгие ограничения на шаг интегрирования по сравнению с одномерным случаем сохраняются. Расщепление второго порядка точности Чтобы обеспечить второй порядок аппроксимации трехмерной схемы, если в качестве базисных были взяты одномерные схемы как минимум второго порядка, используют такую комбинацию F (A1, A2, A3 ) = 1 6 fi (Ai )fj (Aj )fk (Ak ), i=j=k=i (1.29) где суммирование осуществляется по всем шести возможным перестановкам набора цифр (1,2,3). Получившаяся схема наследует максимально допустимый шаг интегрирования от «худшей» базисной схемы fj («худшего» направления j ): < max = max{cmax } j j = c1 max lj j.

Однако расщепление (1.29) обладает рядом недостатков по сравнению с расщеплением (1.27) при практической реализации на компьютере. Можно предложить два алгоритма вычисления (1.29). Первый алгоритм Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ требует вычисления 3! · 3 = 18 действий оператора f и памяти для хранения вспомогательного слоя u : u := f3 un, u := f2 u, un+1 := f1 u. u := f3 un, u := f1 u, un+1 := un+1 + f2 u.... u := f1 un, u := f3 u, un+1 := un+1 + f2 u. u := f1 un, u := f2 u, un+1 := 1 (un+1 + f3 u ). 6 Второй алгоритм позволяет сократить число вычислений f до 15 за счет использования второго вспомогательного слоя u : u := f3 un, u := f2 u, un+1 := f1 u, u := f1 u, un+1 := un+1 + f2 u.... u := f1 un, u := f3 u, un+1 := un+1 + f2 u, u := f2 u, un+1 := 1 (un+1 + f3 u ). 6 Сравните это с тремя вычислениями f в подходе (1.27). Компромиссный вариант Практика показала, что если в (1.29) вычислять только одно любое слагаемое из шести, например, F (A1, A2, A3 ) = f3 (A3 )f2 (A2 )f1 (A1 ), (1.30) то получившаяся схема будет давать результаты очень близкие к способу (1.29). К тому же сохраняется возможность интегрировать с большим шагом по времени, присущим (1.29). Вычисление (1.30): u = f1 un, u = f2 u, un+1 = f3 u, может быть выполнено даже более эффективно, чем вычисление (1.27), за счет того, что разные слои (un, u, u, un+1 ) будут являться лишь псевдонимами одной области памяти.

Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Единственным существенным недостатком такого подхода является несимметрия получающейся схемы, что, впрочем, обычно не критично. Такой способ расщепления использовался, например, в [31], однако, в сочетании с ограничением шага (1.28).

1.6.

Разностные схемы В этом и последующих разделах будем исследовать только одномер ные разностные схемы для системы (1.26), опуская нижний индекс j, имея ввиду построение из них трехмерных схем обозначенными выше способами. Схема Лакса [44] un+1 = un Aun + 2 un, где действие матрицы A рассчитывается по формуле (1.6), 1 un = (un un ), 1 2 +1 1 2 un = (un 2un + un ). 0 1 2 + (1.31) Схема Лакса – Вендроффа [45] Воспользуемся разложением в ряд Тейлора для представления решения на следующем шаге по времени: un+1 = un + Из (1.26) dun 2 d2 un + O( 3 ). + 2 dt 2 dt 2 u d2 u du 2 u =A, = A, dt dt2 2 что после подстановки в ряд Тейлора дает формулу u n+ 1 2 un un A = u A 2 n + O( 3 ), удобную для численного расчета по схеме Лакса – Вендроффа: un+1 = un A(un A2 un ). (1.32) Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Схема Куранта – Изаксона – Риса [46] un+1 = un Aun + 1 ||2 un, (1.33) где || — диагональная матрица, полученная из абсолютных значений. Явный вид действия матрицы из последнего слагаемого нетрудно получить из выражений (1.23) и (1.24): c1 (n · v)n + c2 (n1 · v)n1 + c2 (n2 · v)n2 v 1 = l c1 (T : N00 )N00 + 2c2 ((T : N01 )N01 + (T : N02 )N02 )+. || T +c3 (T : N00 )(N11 + N22 ) Чтобы его упростить и избавиться от вхождения векторов n1, n2, воспользуемся для верхней строки тождеством (1.17) и тождествами (1.15), (1.16), (1.21) для нижней строки. c2 v + (c1 c2 )(n · v)n v 1 || = l (T : N00 )[(c1 2c2 c3 )N00 + c3 I]+. T +c2 [(T · n) n + n (T · n)] Схема Куранта – Изаксона – Риса также может быть записана в форме un+1 = un 1 + un 1 + un, где + un = u n un, +1 0 u n = un un, 0 +, — диагональные матрицы, содержащие только положительные или только отрицательные собственные значения из. Входящие в эту запись схемы комбинации матриц также легко вычислимы, если принять во внимание, что 1 + + 1 = A, 1 + 1 = 1 || :

Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ · n) c2 v + (c1 c2 )(n · v)n + (n · v)I + µ(v n + n v)+ v l 1, = 2 +(T : N00 )[(c1 2c2 c3 )N00 + c3 I]+ T +c2 [(T · n) n + n (T · n)] 1 c2 v + (c1 c2 )(n · v)n (T · n) v l (n · v)I µ(v n + n v)+. 1 + = 2 +(T : N00 )[(c1 2c2 c3 )N00 + c3 I]+ T +c2 [(T · n) n + n (T · n)] Обобщение трех схем Можно предложить следующее однопараметрическое семейство схем, которое включало бы в себя каждую из представленных выше схем: un+1 = un Aun + 1 || 2 un, (1.34) 1 (T где || — результат возведения каждого элемента || в степень. При = 1 получаем схему Куранта – Изаксона – Риса (1.33), при = 2 получаем схему Лакса – Вендроффа (1.32). Случай = 0 порождает схему, максимально похожую на схему Лакса (1.31): || при стремлении к нулю отлична от I лишь для нулевых собственных значениях. В работе [47] была отмечена малая анизотропия решения, полученного по схеме Лакса, в сравнении с другими схемами. У схемы (1.34) при = 0 анизотропия гораздо сильнее. Это обстоятельство натолкнуло на гипотезу, что дополнительное слагаемое 1 (I ||0 )2 un позволяет существенно повысить изотропность схемы, что было подтверждено экспериментально. При нецелых из интервала (0, 2) должны получаться также устойчивые схемы, свойства которых как-то «усредняют» свойства опорных схем, Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ что называется гибридизацией схемы. c2 v + (c1 c2 )(n · v)n v 1 || = l (T : N00 )[(c1 2c2 +2µ c1 )N00 + +2µ c1 I]+. T +c2 [(T · n) n + n (T · n)] (1.35) Для проверки формулы можно убедиться, что подстановка = 2: c2 = 1 + 2µ, c2 = 2 µ, c2 2c2 1 2 c2 = 0, + 2µ дает квадрат (1.6): 2 µ + ( + µ)(n · v)n v v l = A2. 1 ||2 = (T : N00 )I + µ[(T · n) n + n (T · n)] T T Подстановка = 0: v v. 1 ||0 = T:N00 T +2µ [2( + µ)N00 + I] + [(T · n) n + n (T · n)] (1.36) Симметричная схема 4-го порядка точности Будем действовать аналогично выводу схемы Лакса – Вендроффа: dun 2 d2 un 3 d3 un 4 d4 un + + + + O( 5 );

2 3 4 dt 2 dt 6 dt 24 dt 2 2 3 3 4 u du du d4 u du 2 u 3 u 4 u = A, =A, = A, =A ;

dt dt2 2 dt3 3 dt4 4 1 4 un un 1 3 un 1 2 un n + O( 3 );

= u A A A A 2 3 4 2 6 24 un+1 = un + un+1 = un A(un A(2 un A(3 un A4 un ))). (1.37) un+ Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Здесь разности аппроксимируют соответствующие производные с точностью 4-го порядка: 1 (2un + 16un 16un + 2un ), +2 +1 1 2 24 1 2 un = (un + 16un 30un + 16un un ), +2 +1 0 1 2 24 1 3 un = (2un 4un + 4un 2un ), +2 +1 1 2 24 1 4 un = (un 4un + 6un 4un + un ). +1 0 1 2 24 +2 un = 1.7.

Сеточно-характеристические схемы Среди всех разностных схем для гиперболических уравнений (1.26), т.е. таких, матрица которых имеет полный набор собственных векторов и только действительные собственные значения: u = 0, u + 1 (1.38) принято выделять так называемые сеточно-характеристические схемы [22], которые существенно учитывают особенности таких уравнений. Специфика гиперболических уравнений состоит в характеристических гиперповерхностях (при одной пространственной координате — прямых), на которых уравнения не содержат выводящих производных. Видно, что система (1.38) при замене переменных на самом деле распадается на ряд независимых скалярных уравнений переноса v = 0. (v u) v+ Уравнения переноса решаются весьма просто (рис. 1.1). Из того узла m временного слоя n + 1, в котором требуется получить решение, опускаются характеристики. Из точки пересечения характеристики со слоем n n+1 соответствующая компонента вектора v переносится в точку m : n+1 n vi (m ) = vi (m i ).

Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ n+1 t l>0 m-1 lt Рис. 1.1. Сеточно-характеристический принцип расчета узла.

n+1 n В последующем мы будем опускать m, записывая просто vi = vi (i ).

l=0 m l<0 m+1 n Если характеристика не попадает точно в расчетный узел, то применяются те или иные способы реконструкции решения в данной точке (интерполяция). Выбор способа реконструкции решения на предыдущем слое обуславливает степень аппроксимации сеточно-характеристической разностной схемы. В работе [22] приводится вывод сеточно-характеристических схем первого (Куранта – Изаксона – Риса) и второго (Лакса – Вендорффа) порядков. После того, как все компоненты v перенесены, восстанавливается само решение: un+1 = 1 v n+1. В некоторых программных реализациях сеточно-характеристического метода непосредственно производится обращение матрицы или эквивалентное ей решение системы линейных уравнений. Здесь же будет предложен другой подход. Подставим одно из последних выражений в другое ··· ··· n+1 1 n 1 i n = vi (i ) = · u (i ) = i · un (i ),i = u i ··· ··· = i,i i un (i ), где i — i-ая строка,,i — i-ый столбец 1. Обозначив Xi =,i i, Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ окончательно получаем un+1 = i Xi un (i ).

(1.39) Альтернативно можно собрать вместе сумму всех элементов, отвечающих одному и тому же собственному значению: Xc = i =c,i i.

Поскольку I = 1 то Xi = I i =0 i = T = i T Xi, Xi, и можно придать (1.39) форму часто более удобную для применения: un+1 = un + i = Xi [un (i ) un ].

(1.40) Еще раз подчеркнем, что полученная формула справедлива для любой сеточно-характеристической схемы, не зависимо от того, какого порядка реконструкция функции на предыдущем слое по времени будет использоваться. Характеристические матрицы для системы уравнений упругости Вывести действие матриц X на вектор не составит труда, если обратится к выражениям (1.25) и (1.23): v n 1 1 = N00 : T) X±c1 = (n · v 2 c1 T [(c1 c3 )N00 + c3 I] 1 (n · v)n c1 (N00 : T)n 1, = T:N00 2 (n · v)[(c1 c3 )N00 + c3 I] + (I + 2µN00 ) +2µ Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ v n1 1 1 + N01 : T) X±c2 = (n1 · v 2 c2 2c2 N01 T n2 1 1. N02 : T) + (n2 · v 2 c2 2c2 N02 Чтобы во второй формуле оставить только n и N00, воспользуемся тождествами (1.17), (1.18), (1.20), (1.21). Тогда 1 v (n · v)n c2 [T · n (N00 : T)n] v 1 ±c2 X =. c2 [n v + v n 2(n · v)N00 ]+ 2 T +(T · n) n + n (T · n) 2(T : N00 )N00 Сумма характеристических матриц для положительных (отрицательных) собственных значений определяется суммированием двух последних выражений: v v X± = (X±c1 + X±c2 ) = T T 1 11 1 (1.41) v c2 T · n ± ( c2 c1 )(N00 : T)n 1 = (n · v)[(c1 2c2 c3 )N00 + c3 I] c2 (n v + v n)+. 2 T:N00 + +2µ [I 2( + µ)N00 ] + (T · n) n + n (T · n) Нетрудно проверить, что X + X+ = 1 ||0, см. (1.36).

1.8.

Расчет на границе области интегрирования Формулы (1.39) и (1.40) не могут быть непосредственно применены вблизи границы, поскольку в них входят величины un (i ), которые находятся за пределами области интегрирования. Корректное решение задачи на границе требует постановки граничных условий в количестве равном числу выходящих характеристик [22]. Для левой границы области выходящими являются характеристики, отвечающие положительным собственным значениям, для правой границы — отрицательным (рис. 1.1). Общий Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ вид линейных граничных условий имеет вид Bun+1 = b, где B — прямоугольная матрица, b — вектор правых частей. Разделим все множество собственных значений A на две группы по признаку попадания характеристик внутрь области на предыдущем временном слое: int и out. Будем считать, что выражение (1.39) сохраняет справедливость и для граничных узлов, если величины un (i ), i out специально подобраны так, что выполнено (1.42): B i (1.42) Xi un (i ) = b, что эквивалентно Xi un (i ) + i out i in B i,i = b, где i — некоторые коэффициенты, при которых справедлив знак равенства. Обозначим первую сумму как uin, а все столбцы,i и коэффициенты i из второй суммы соберем в прямоугольную матрицу,out и вектор соответственно: B uin +,out = b. Значения вектора вычислять нет необходимости, поэтому исключим его из последнего выражения: = B,out (b Buin ).

Подставляя его обратно в (1.39), находим un+1 = uin +,out B,out (b Buin ).

(1.43) В итоге получилась обобщенная запись метода расчета приграничного узла (хоть одна выходящая характеристика из которого не попадает внутрь Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ области интегрирования на предыдущем шаге по времени). Формула справедлива для любого порядка точности реконструкции функции на слое n и для любого набора граничных условий {B, b}. Выражение (1.43) помимо прочего задает критерий корректности постановки того или иного граничного условия, а именно: невырожденность матрицы B,out. В последующих подразделах будет показано, как эта запись упрощается для уравнений линейной упругости с определяющей матрицей (1.6) и конкретного вида граничных условий. 1.8.1. Заданная внешняя сила Данное условие заключается в том, что на границе тела задана поверхностная плотность силы воздействия со стороны внешних тел. В частном случае эта сила может быть равна нулю, тогда говорят о свободной границе тела. Условие формализуется следующей записью: Tp = f, где p определяет перпендикуляр (единичной длины) к поверхности тела, направленный наружу, f — поверхностная плотность внешних сил. Данная запись содержит три скалярных условия, которые компенсируют отсутствие данных для трех выходящих характеристик. Приводя граничные условия к виду (1.42), получаем v B = Tp, T Что после подстановки в (1.43) дает un+1 = uin,out B,out где uin = {v in, Tin }.

b = f.

(Tin p f ), (1.44) Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Действие,out устанавливается из (1.24): w 1 1 w 1 n + w2 n 1 + w 3 n 2,,out w2 = 2 [w1 ((c1 c3 )N00 + c3 I) + 2c2 (w2 N01 + w3 N02 )] w3 w1 B,out w2 = [w1 ((c1 c3 )N00 + c3 I) + 2c2 (w2 N01 + w3 N02 )]p, 2 w3 где знак минус относится к левой, а плюс — к правой границе области. Для раскрытия (1.44) необходимо в выражении 2 w = B,out z 2 B,out w = z 2 вычислять w для произвольного z = Tin p f, коэффициент выбран 2 для упрощения последующих формул. Помножим z = B,out w после довательно на n, n1, n2 : n · z = w1 c1 (n · p) + w2 c2 (n1 · p) + w3 c2 (n2 · p), n1 · z = w1 c3 (n1 · p) + w2 c2 (n · p), n2 · z = w1 c3 (n2 · p) + w3 c2 (n · p). В результате w1 = (n · p)(n · z) (n1 · p)(n1 · z) (n2 · p)(n2 · z), c1 (n · p)2 c3 [(n1 · p)2 + (n2 · p)2 ] w2 = n1 · b, w3 = n2 · b, 1 b(z) (z w1 c3 p). (1.45) c2 (n · p) Практическая реализация в программе иногда требует, чтобы формулы не содержали векторов n1 и n2, в выборе которых есть некоторая степень свободы. Чтобы исключить их, распишем скалярное произведение произвольных векторов a и b в ортогональной системе координат с базисом (n(0), n1, n2 ): (ni · a)(ni · b) = (a · b), i Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ что дает w1 (z) = 2(n · p)(n · z) (p · z). (c1 + c3 )(n · p)2 c3 p (1.46) Окончательно, учитывая тождество (1.18) для вектора b, 1 [(w1 b · n)n + b] 1,,out B,out z= [(c1 c3 )w1 2c2 (n · b)]N00 + c3 w1 I + c2 (b n + n b) где b и w1 есть функции z = Tin p f : (1.45) и (1.46) соответственно.

Поверхность тела лежит в координатной плоскости сетки Дальнейшие упрощения возникают, когда две прочие линии сетки по отношению к той, для которой строятся расчетные формулы, точно совпадают с поверхностью тела. Это имеет место, например, при моделировании тела лагранжевой сеткой, узлы которой сдвигаются вместе с границей. Тогда нормаль к границе будет совпадать с перпендикуляром к двум прочим линиям сетки, то есть с вектором биортогонального базиса. Без ограничения общности можно считать, что p = n, где знак минус относится к левой, а плюс — к правой границе области. Тогда, обозначив z = Tin n f, 1 w1 = (n · z), c1 b= 1 c3 (n · z) n, z c2 c1 1 c3 n · b = (n · z) 1. c2 c Окончательно (1.44) приобретает такой покомпонентный вид для границы, перпендикуляр которой коллинеарен вектору биортогонального базиса n, v n+1 = v in T n+1 in z·n (I 2( + µ)N00 ), = T ± [z n + n z] ± + 2µ 1 1 z+ c 1 1 c2 c (z · n)n, (z = Tin n f ) (1.47) где, напомним, верхний знак относится к левой, а нижний — к правой границе области интегрирования.

Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ 1.8.2.

Заданная скорость границы Данное условие заключается в том, что на границе тела задана скорость движения тела, например, это может быть скорость внешнего тела, контактирующего с данным. В частном случае скорость может быть равна нулю, тогда граница тела будет неподвижна. Условие формализуется следующей записью: v =V, где V — заданная скорость. Данная запись содержит три скалярных условия, которые компенсируют отсутствие данных для трех выходящих характеристик. Приводя граничные условия к виду (1.42), получаем v B = v, b = V. T Действуя аналогично с условием заданной силы на границе, последовательно находим 2B,out w = w1 n + w2 n1 + w3 n2, обозначаем «невязку» z = v in V, разрешаем n·z w1 1 1 z = n1 · z, B,out w2 = 2 w3 n2 · z подставляем,out B,out n·z 1 z = 2,out n1 · z = n2 · z z. = [(n · z)((c1 c3 )N00 + c3 I) + 2c2 ((n1 · z)N01 + (n2 · z)N02 )] Окончательный покомпонентный вид корректировки на границе с за Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ данной скоростью, учитывая тождество (1.18) для вектора z, имеет вид v n+1 = V, Tn+1 = Tin ± [(n · z)((c1 2c2 c3 )N00 + c3 I) + c2 (z n + n z)], (1.48) где z = v in V и, напомним, верхний знак относится к левой, а нижний — к правой границе области интегрирования. 1.8.3. Смешанные условия Первое смешанное условие заключается в том, что заданы нормальная скорость движения границы и тангенциальная составляющая силы. Такое условие возникает в случае контакта с другим телом. Тангенциальная составляющая силы равна нулю в случае отсутствия трения при контакте и пропорциональна скорости скольжения в противном случае. Обозначим заданную касательную силу символом f, компонента этого вектора в направлении нормали может быть произвольной, заданную нормальную скорость — Vp. Тогда смешанное условие формализуется записью Tp = f, v · p = Vp, f = f (f · p)p + (f · p)p, где (f · p) предстоит установить. Умножим первое уравнение в (1.47) на n: (v n+1 · n) = (v in · n) 1 (z · n), c1 (z = Tin n f ) (1.49) левая часть равенства есть Vp, следовательно (f · n) = Tin : N00 c1 [(v in · n) ± Vp ], откуда вычисляется поверхностная плотность силы, действующей в данной точке f = f + n[c1 Vp + (Tin · n c1 v in f ) · n], (1.50) Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ верхний знак относится к левой, а нижний — к правой границе области интегрирования. Расчет границы, на которой выставлено первое смешанное условие, ведется по формуле (1.47), в которую подставляется f, представленное в последней формуле. Для полноты картины опишем другое смешанное условие, когда заданы тангенциальная скорость движения границы и нормальная составляющая силы. Если обе заданные величины равны нулю, мы имеем дело с «шершавой» свободной границей, где точке границы разрешено двигаться только перпендикулярно, но зато внешняя нормальная сила равна нулю. На практике потребность во втором смешанном условии чаще связана с задачами, описанными в подразделе 1.8.4. Обозначим заданную касательную скорость символом V (компонента этого вектора в направлении нормали может быть произвольной), нормальное давление на границе — fp. Тогда второе смешанное условие формализуется записью p · Tp = fp, v = V (V · p)p + (v · p)p, где (v · p) устанавливается из (1.49): v = V + n (v in V ) · n ± 1 (n · Tin n fp ), c1 (1.51) верхний знак относится к левой, а нижний — к правой границе области интегрирования. Расчет границы, на которой выставлено второе смешанное условие, ведется по формуле (1.48), в которую подставляется v, представленная в последней формуле. 1.8.4. Условия поглощения и симметрии Рассмотрим пример моделирования распространения плоских волн от дневной поверхности (геол.) вглубь некоторой породы, в котором нас Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Рис. 1.2. Пример, когда некоторые границы области интегрирования (пунктир) не являются физическими границами моделируемой сплошной среды.

интересует отражение от неоднородностей среды, также наблюдаемое на поверхности (рис. 1.2). Область моделирования показана пунктирной линией. При этом верхняя ее граница совпадает с дневной поверхностью, и там логично поставить условия нулевой внешней силы. Тогда как боковые и нижняя границы проведены произвольно, но для численного расчета на них также требуется постановка граничного условия. Для того чтобы боковые границы не вносили искажений в результат наблюдения, они должны выбираться по возможности сильно удаленными от неоднородностей среды, являющихся источником отраженных сигналов. В этих условиях состояние среды на этих границах будет зависеть только от возмущения, исходящего от дневной поверхности. Поскольку это возмущение имеет плоский фронт, то картина слева и справа от вертикальной границы будет симметричной. Граничные условия на боковых поверхностях существенно зависят от вида падающей волны: продольная она или поперечная. Для продольной волны составляющая скорости в направлении нормали к данной границе должна быть нулевой. Из-за того что главные оси тензора напряжений должны совпадать с боковой границей, тангенциальная составляющая силы, действующей между частями среды, условно разделенными данной границей, также обращается в нуль. Поэтому здесь мы можем производить расчет по правилам смешанной границы (подраздел 1.8.3) с условиями f = 0, Vp = 0.

Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Для поперечной волны, наоборот, должны выполняться равенства fp = 0, V = 0.

Эти условия становятся очевидными, если обратиться к виду собственных векторов (1.25), одному из которых (в зависимости от типа и направления движения волны) пропорциональны скорости и напряжения во всех точках падающей плоской волны (вектор p направлен вдоль волнового фронта). Для нижней границы ситуация иная. Здесь требуется, чтобы граница не отражала падающий на нее сигнал, а полностью поглощала его, поскольку предполагается, что реальная граница тела, от которой возможно отражение, находится на бесконечности. За волну, распространяющуюся от границы тела, отвечают компоненты i u состояния граничного узла u, относящиеся к выходящим характеристикам. Поэтому условие на поглощающей границе, приведенное к виду (1.42), примет вид out u = 0, где out — прямоугольная матрица, состоящая из строк (1.23), относящихся к выходящим характеристикам для данной границы. Однако реализовать полное поглощение энергии таким способом возможно лишь для волн распространяющихся в направлении p. Для прочих волн гарантируется лишь частичное поглощение, поэтому условия поглощения желательно комбинировать с удалением такой границы на возможно большее расстояние. Второй этап расчета для поглощающей границы определяется из выражения (1.43): un+1 = (I,out out )uin. Произведение матриц,out out для p = n равняется X± (1.41). Отсюда T Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ покомпонентная запись (1.43): 1 in 1 in 11 1 v n+1 = v± T ·n (N00 : Tin )n 2 c2 c2 c1 1 Tn+1 = Tin + {±(n · v in )[(c1 2c2 c3 )N00 + c3 I] ± c2 (n v in + v in n) 2 in T : N00 [I 2( + µ)N00 ] [(Tin · n) n + n (Tin · n)]}. + 2µ (1.52) Если сравнить эти выражения с (1.47) и (1.48), то окажется, что нет необходимости реализовывать (1.52) в программе отдельно. Достаточно взять полусумму результатов, полученных указанными корректорами для свободной (f = 0) и фиксированной (V = 0) границ. После достижения времени, необходимого излученной с дневной поверхности плоской волне полностью пройти нижнюю границу области моделирования, условия на боковых границах желательно изменить на условия поглощения. 1.8.5. Решение на границе при наличии правой части В этом подразделе рассматривается расчет граничного узла для неоднородной системы уравнений вида (1.26), которая содержит правую часть r, являющуюся линейной функцией решения: r = Ru. Обычно в таком случае прибегают к методу расщепления. Переход от un к un+1 осуществляется в два этапа через вспомогательный слой u : un u un+1. В качестве одного из этапов выступает интегрирование однородной системы уравнений (1.26), решение которого уже было выписано (1.43). Другим этапом является интегрирование неоднородного обыкновенного дифференциального уравнения: u = r.

Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Решение этого уравнения можно представить как действие матрицы F, зависящей от шага интегрирования и R: u(t + ) = Fu(t). Например, для неявного метода первого порядка F = (I R)1. Если на первом этапе (un u ) решается однородное уравнение, а на втором (u un+1 ) — учитывается наличие правой части, то un+1 = Fuin + F,out B,out (b Buin ).

(1.53) Если поменять порядок следования этапов: сначала «добавляем» правую часть, затем решаем однородное уравнение, то un+1 = Fuin +,out B,out (b BFuin ).

(1.54) В качестве третьего способа рассмотрим более сложную конструкцию. Пусть интегрирование правой части уравнения составляет второй этап, но при этом потребуем выполнения граничных условий для un+1 (1.42). Повторяя действия, изложенные в ходе получения (1.43), запишем систему алгебраических уравнений u = uin +,out, un+1 = Fu, Bun+1 = b. После исключения из нее u и : un+1 = Fuin + F,out BF,out (b BFuin ).

(1.55) Выражения (1.53), (1.54), (1.55) стремятся друг к другу при F I. Результаты их применения также оказываются очень похожими. Однако расчеты контактных границ при осевой симметрии задачи [48], проводимые в цилиндрической системе координат (переход в которую приводит к появлению правых частей в уравнениях), показали, что только (1.55) полностью исключает проявление неустойчивости вблизи контактных поверхностей.

Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ 1.9.

Контакт между двумя телами va, Ta vb, Tb p a из областей.

b Рис. 1.3. Контакт двух областей. Расчетный узел содержит по набору переменных для каждой Предположим, что два тела a и b, возможно из различных материалов, соприкасаются в рассматриваемой точке вдоль какой-то поверхности с нормалью p, внешней по отношению к телу a. Кроме того, два граничных узла из сеток, построенных в телах a и b, полностью совпадают (рис. 1.3). Необходимо для данных узлов построить корректоры в стиле раздела 1.8, которые бы компенсировали недостаток характеристик, попадающих в свое тело, и учитывали особенности того или иного типа контакта между телами. Далее рассматриваются различные типы контактов. 1.9.1. Полное слипание При полном слипании двух тел, скорости соприкасающихся граничных точек должны совпадать, а внешние силы, действующие на каждое тело, равны с точностью до инвертирования направления: va = v b = V, fa = fb. (1.56) Внешние силы, действующие на тела на следующем шаге, равны соответственно fa = Tn+1 p, a fb = Tn+1 (p). b Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Домножим второе выражение в (1.48) скалярно на p: fa = Tn+1 p = Tin p a {(p · za )(ca1 ca2 )p + ca2 za }, a a in (za = va V ) in fb = Tn+1 p = Tin p + b {(p · zb )(cb1 cb2 )p + cb2 zb }. (zb = vb V ) b b (1.57) Расписывая разность (Tn+1 Tn+1 )p = 0, a b не составляет большого труда определить общую скорость обоих контактирующих узлов: V= 1 in in {a [(p · va )(ca1 ca2 )p + ca2 va ]+ a ca2 + b cb in in + b [(p · vb )(cb1 cb2 )p + cb2 vb ] (Tin Tin )p a b a (ca1 ca2 ) + b (cb1 cb2 ) in in (p · [a ca1 va + b cb1 vb (Tin Tin )p])p}. a b a ca1 + b cb1 (1.58) Как и следовало ожидать, получившееся выражение инвариантно относительно переименования областей a b c одновременным изменением знака у p. После того как скорость найдена, значения в каждом из узлов корректируются независимо по формулам (1.48). 1.9.2. Свободное скольжение При свободном скольжении двух тел вдоль общей поверхности раздела должны совпадать лишь нормальные составляющие скорости, и нормальные компоненты силы должны компенсировать друг друга, в то время как тангенциальные компоненты силы обращаются в нуль: va · p = vb · p = Vp, a b fp = fp, fa = fb = 0.

(1.59) Внешние силы, действующие на каждое из тел, определяются из (1.50):

a in fp = a ca1 Vp + (Tin p a ca1 va ) · p a b in fp = b cb1 Vp + (Tin p b cb1 vb ) · p. b (1.60) Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Из того, что сумма этих сил должна равняться нулю, следует Vp = 1 in in a ca1 va + b cb1 vb (Tin Tin )p · p. a b a ca1 + b cb1 (1.61) Заметим, что это выражение совпадает с нормальной составляющей скорости при полном слипании (1.58). Таким образом, нормальная скорость движения контакта, а следовательно, и нормальная составляющая силы между телами не зависит от типа контакта. Если свойства упругих материалов обоих тел совпадают, то Vp = 1 in 1 in va + vb (Tin Tin )p · p. b 2 c1 a Найденная нормальная скорость движения контактной границы может быть подставлена в выражение для расчета границы при поставленном смешанном условии (1.50). Так же, как и для случая полного слипания, коррекция контактирующих узлов далее ведется независимо.

1.10.

Интегрирование уравнений акустики Важным частным случаем уравнений упругости (1.2) являются уравнения акустики. Они возникают, когда тензор напряжений имеет шаровую, «гидростатическую» форму: T = pI и, следовательно, µ = 0. Упрощенная система (1.2) получается v = p, p = ( · v). (1.62) Приведем далее вид наиболее важных формул, полученных в результате исследования уравнений упругости и численных методов их интегрирования. Выражение (1.4) остается справедливым, если под вектором u понимать u = {v, p} = {v1, v2, v3, p}, T T Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ а действие матрицы Aj на произвольный вектор p nj v. Aj = lj (nj · v) p Комбинации констант: c1 = c3 =, c2 = 0.

У матриц A два нулевых и два ненулевых собственных значения: = l diag{c1, c1, 0, 0}, p n · v + c1 w1 w 2 n · v cp v 1 = =, w 3 n1 · v p w4 n2 · v v (w + w2 )n + w3 n1 + w4 n2 = 1 w = 1 1. 2 c1 (w1 w2 ) p Важная комбинация матриц (n · v)n v 1 || = c 1 p p позволяет воспроизвести запись любой из приведенных разностных схем.

1.11.

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

Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Система дифференциальных уравнений (1.2), а также запись (1.6) остаются в силе, при условии того, что вектор скорости и тензор напряжений теперь включают меньшее количество компонент: u = {v, T} = {v1, v2, t11, t12, t22 }. У матриц A две пары ненулевых и одно нулевое собственное значение: = l diag{c1, c1, c2, c2, 0}, 1 w N :T 1, w2 n · v c1 00 v 1 w3, w4 = = n1 · v c2 N01 : T, T w5 (N11 +2µ N00 ) : T (w1 + w2 )n + (w3 + w4 )n1 v 1 = 1 w = (w2 w1 )[(c1 c3 )N00 + c3 I]+, 2 T +2c2 (w4 w3 )N01 + 2w5 (I N00 ) 0 n n1 1. 1 = 2 [(c1 c3 )N00 + c3 I] 2c2 N01 2(I N00 ) Здесь для упрощения записи элементы (строки и столбцы) матриц, отличающиеся лишь знаком, объединены, причем верхний знак в соответствует первому элементу, а нижний — последующему. Запись комбинации матриц 1 || (1.35) при этом также сохраняется.

T T 1.12.

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

Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ Для расчета во внутренних неподвижных узлах необходимо учитывать наличие конвективных членов, появляющихся из-за движения тела относительно узлов сетки: u x u u + = + v · u. u= t t x t Тогда, по-прежнему, будет справедлива запись системы уравнений в форме (1.4), если только переопределить входящие в нее матрицы: Aj = lj (nj · v)I + Aj, где lj (nj · v) — контравариантные координаты скорости тела в данной точке относительно выбранного базиса {j } и взаимного к нему {lj nj }, Aj определяются выражением (1.6). Максимальные абсолютные величины собственных значений Aj равны cj = c1 + |lj (nj · v)|. 1 Именно эти величины будут определять максимально допустимый шаг интегрирования.

x1 b+ x1 b x2 x2 i n2 n1 O b Рис. 1.4. Расчет граничного узла при совместном использовании подвижных маркеров (крестики) с регулярной квадратной эйлеровой решеткой (кружки). Значения во внутренних узлах, отмеченных полыми кружками, определяются интерполяцией.

Алгоритмы расчета неподвижного внутреннего узла и подвижного граничного узла не отличаются от изложенных в предшествующих разделах. Для внутренних узлов следует иметь ввиду, что они могут располагаться сколь угодно близко к линии границы. Поэтому непосредственный Глава 1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ УПРУГОСТИ расчет таких узлов будет диктовать стремящийся к нулю шаг интегрирования. Из-за чего внутренние узлы имеет смысл разделить на два класса: 1) отстоящие от границы дальше, чем от соседей по решетке и 2) прочие приграничные внутренние узлы (рис. 1.4). Для первых — напрямую применять один из разностных методов, а для вторых — получать значения интерполяцией по ближайшим уже рассчитанным внутренним и граничным узлам. Здесь также имеет смысл рассмотреть каким образом можно аппроксимировать частные производные при совместном использовании подвижных маркеров с регулярной квадратной эйлеровой решеткой. Допустим, необходимо получить значения, отвечающие следующему моменту времени для граничного узла b (рис. 1.4). Если введено такое направление обхода границы, что внутренность тела остается слева, то можно обозначить предшествующий граничный узел как b, а последующий — как b+. Также для расчета потребуются данные одного из внутренних узлов тела, который мы обозначим i. Перейдем в систему координат (1, 2 ), которая определяется своими единичными векторами (в мировой системе координат их длина не единичrb+ rb, 2 = ri rb, 2 где символом r обозначено положение узлов. Нуль O в этой системе коорди1 = нат совпадает с точкой пересечения двух прямых: одной проходящей через узлы b, b+, другой — b, i. Предполагается, что мелкость сетки выбрана так, что узел b, точка O и центр отрезка [b, b+] практически совпадают, то есть можно пренебречь изломом границы около узла b и считать плотность узлов на границе равномерной. Вдоль линий новой системы координат легко определить первые и вторые разности, которые подставляются в формулы разностных схем (см. раздел 1.6): ub+ ub ub+ 2ub + ub, 21 ub =, 1 ub = 2 2 ная) + u b = u i ub. Глава 2 Построение нерегулярной треугольной сетки При численном моделировании для описания состояния тела обычно вводится некоторый конечный набор точек в теле, в которых и только в них запоминаются значения переменных. В остальных точках решение можно приближенно восстановить из значений, хранимых в нескольких близлежащих точках набора, откуда следует важность отыскания необходимого количества ближайших точек. Под сеткой понимают комбинацию из введенного набора точек (далее узлы сетки) и дополнительной информации о их соседстве. Численный метод расчета ответственен за вычисление значений в узлах сетки в следующий момент времени по их прежним значениям. В вычислительной математике принято разделять сетки на структурированные (регулярные) и неструктурированные (нерегулярные). В структурированных сетках информация о соседстве существенно более проста, примером может служить прямоугольная сетка, использованная во множестве работ, где для каждого внутреннего узла можно узнать ссылки на четырех соседей вдоль линий решетки. В неструктурированных сетках ограничения на количество соседей и положения узлов существенно слабее, зато информация о соседстве устроена значительно сложнее. Потребность в моделирование тел со сложной геометрией, которая задана изначально, либо формируется в результате деформаций в процессе Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ расчета, возникает все чаше. Однако даже для таких задач исследователи зачастую продолжают использовать регулярные сетки. Существенным недостатком является возможность их использования лишь для ограниченных классов геометрий тел, поскольку форма границ тела и наличие в нем полостей вводят ограничение на возможное положение узлов сетки. Как правило, для автоматического построения структурированной сетки требуется односвязность области. Кроме того, даже если регулярную сетку удается построить, бывает невозможно гарантировать одинаковую плотность узлов во всех частях сетки, а также малую кривизну сеточных линий. Например, в круговой области деформированная решетка сгущается и искривляется при приближении к четырем точкам границы [49]. Эти особенности приводят к деградации точности решения для численных методов, основанных на использовании решеток. Кроме того, изменение оптимального положения узлов решетки при деформации тела зачастую требует постоянного решения ресурсоемкой многомерной задачи минимизации некоторого функционала [35]. Также из регулярной решетки невозможно произвольно выкидывать отдельные узлы, скажем, для уменьшения их плотности в местах сильного сжатия и поддержания стабильного шага интегрирования. Фрагментация тел возможна лишь вдоль линий или плоскостей решетки. В работе [17] более подробно обсуждаются проблемы, связанные с использованием структурированных сеток, и описываются различные базовые подходы решения динамических задач на неструктурированных сетках. Данная глава посвящена алгоритмам работы с треугольной неструктурированной сеткой. Описывается формат представления сетки в памяти компьютера — двусвязный список ребер, обеспечивающий константное время выполнения большинства элементарных операций поиска и модификации сетки. Приводятся подробные изложения алгоритмов, необходимых для создания начальной сетки исходя только из контуров тела, а также Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ методы ее перестройки, обеспечивающие требуемые свойства сетки при деформировании тела в процессе моделирования. Построенная сетка является триангуляцией Делоне. Доказывается ряд утверждений касательно размеров ячеек получаемой в результате работы алгоритмов сетки. Формулируется ограничение на максимальный шаг интегрирования, при котором гарантируется невырожденность сетки. Глава завершается несколькими примерами сеток, полученных в результате выполнения описанных алгоритмов.

2.1.

Представление триангуляции в программе Прежде чем приступить к описанию алгоритмов построения и вы равнивания сетки необходимо объяснить, в виде каких структур данных триангуляция тела хранилась в программе. Выбор формата данных определяет временные характеристики и затраты памяти алгоритмов, обрабатывающих сетку, поэтому чрезвычайно важно выбрать такое представление сетки в памяти, чтобы максимально облегчить последующие операции с ней. 2.1.1. Наиболее компактный формат Часто для описания тела, разбитого на треугольники, используют два массива: один, ставящий в соответствие номеру узла его координаты, и другой, формирующий треугольники, используя номера трех его вершин. Элементы этих массивов на языке Си можно объявить следующим образом: struct node_data { float x, y;

//position... //other node’s data };

struct tri_data { Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ unsigned int nodes[3];

//vertices... //other triangle’s data };

Подобное описание обладает рядом положительных качеств: • Минимальность и простота. Трудно придумать другое описание, которое бы занимало существенно меньше места и было бы столь же очевидным в интерпретации. Ввиду этого оно идеально подходит для долговременного хранения треугольных сеток. • Удобство. Ряду алгоритмов машинной графики необходимо лишь быстро перебирать треугольники и в каждом из них учитывать только значения локальных вершин. К этому ряду принадлежит, например, алгоритм построения изолиний marching cubes. • Широкое распространение. В таком виде данные принимают многие приложения, генерирующие сетку или использующие ее. Однако бывает, что выполнение некоторых операций при таких структурах данных затруднено необходимостью затрачивать долгое время (пропорциональное числу узлов или треугольников) на поиск: • Модификация сетки: введение или исключение новой вершины с инцидентными ей ребрами и смежными треугольниками, переключение ребра внутри пары треугольников и т. д. • Однократное рисование каждого ребра сетки. Если выводить ребра всех треугольников, то каждое ребро за исключением граничных будет выведено дважды (по количеству треугольников, которые ее разделяют). • Переход от одного треугольника к соседнему через общее ребро. Это может потребоваться, например, для построения гладких изолиний, когда нужно усреднять нормали к линиям по соседним треугольникам.

Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ • При реализации численных методов и переинтерполяции на новую сетку важен поиск треугольника, в который попадает точка, находящаяся от некоторого узла на определенном смещении. Ввиду этих потребностей после загрузки сетки с долговременного носителя в память программы, ее внутреннее представление должно быть расширено, но желательно все же заботиться об экономии памяти и не хранить те данные, которые могут быть вычислены за константное время. 2.1.2. Расширенные структуры данных для ускорения поиска Поскольку треугольная сетка является частным случаем планарного графа, то в качестве основы ее представления в данной работе была выбрана структура данных, предназначенная для программного описания подразделений плоскости, наиболее близким известным аналогом которой является двусвязный список ребер (doubly-connected edge list) [42], в который были внесены некоторые модификации, представленные здесь. Несмотря на наличие в названии слова «список», сетка представлялась главным образом с помощью нескольких массивов, а не списков, которые содержат объекты сетки без пропусков для рационального расходования памяти. Ссылкой на объект сетки в таком случае является просто индекс в соответствующем массиве. Хранимыми объектами являются узлы, полу-ребра, ребра, треугольники. Элементы этих массивов можно объявить следующим образом: struct node_data { float x, y;

//position int he;

//index of an incident half-edge //or -1 if the node has no incident edges... //other node’s data };

Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ struct half_edge_data { int tri;

// index of triangle to the left //or -1 if no left triangle exists unsigned int node;

// index of base node unsigned int next;

// index of next half-edge (of the same // node) in counter-clockwise order... //other half-edge’s data };

struct edge_data {... //other edge’s data };

struct tri_data { unsigned int nodes[3];

// indices of vertices... //other tri’s data };

Из всех объектов наиболее сложным образом представлены грани: записями сразу в двух массивах. Каждое ребро условно разбивается на две половины — полу-ребра, которые в отличие от целого ребра являются направленными, и для них определены понятия «слева» и «справа». Каждое полу-ребро хранит ссылки на смежный треугольник, находящийся слева, базовый узел, из которого полу-ребро «растет» и следующее полу-ребро данного узла в направлении против часовой стрелки. Ссылка на парное (twin) полу-ребро не хранится, а вычисляется, поскольку они всегда возникают и исчезают парами, и в массиве всегда располагаются в смежных позициях. Количество элементов в массиве ребер вдвое меньше, чем в массиве полу-ребер, а индексы ребра e и его частей h0, h1 связаны простым соотношением: h1 1 = h0 = 2e. Ребро не хранит никакую информацию, необходимую для поиска в сетке, а необходимо лишь как контейнер для ассоциированных с ребром данных.

Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ Каждый узел помимо своих координат содержит ссылку на одно из инцидентных полу-ребер. Треугольник, по-прежнему, представляется ссылками на три свои вершины. Однако в результате экспериментов с такой структурой данных выяснилось, что все же удобнее вместо ссылок на вершины хранить ссылки на полу-ребра, оставляющие треугольник слева, или даже всего на одно из таких полу-ребер для экономии памяти. Помимо массивов также поддерживался набор индексов граничных полу-ребер, представляющих каждое одну из границ в сетке. Это позволяет осуществлять перебор только граничных вершин, ребер или треугольников. Хотя необходимость поддержания этого набора может замедлить процедуры модификации сетки. Отдельного упоминания заслуживает процедура удаления объектов из сетки. Поскольку пропуски в массивах не допускаются, то удаленный элемент заменялся последним элементом массива с сокращением на единицу размера последнего. Поэтому удаление элементов приводит к недействительности внешних ссылок на оставшиеся элементы этого массива. В большинстве случаев алгоритмы могут быть легко модифицированы с учетом этого обстоятельства. Например, процедура сканирования всех элементов массива в порядке возрастания индексов с удалением объектов, удовлетворяющих некоторому критерию, должна быть видоизменена так, чтобы после удаления повторно проверялся элемент с тем же индексом. Однако иногда бывает иначе, когда потеря состоятельности ссылок недопустима. Для таких случаев специально были разработаны процедуры удаления с отложенной перенумерацией оставшихся элементов.

Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ 2.2.

Триангуляция невыпуклого многоугольника с полостями Метод построения сетки, приведенный далее, основывается на име ющейся триангуляции тела. В самом начале расчета программа обладает информацией лишь о контурах тела, на основании которой ей необходимо породить какую-либо триангуляцию области. В самом общем случае форма тела задается ломанными контурами своих единственной внешней и возможно нескольких внутренних границ. Таким образом, появляется потребность триангулировать произвольный невыпуклый многоугольник, содержащий множество полостей. Количество полостей в проведенных расчетах достигало нескольких тысяч, поэтому этот аспект подчеркивается не случайно. Можно пойти по простому пути, использованному в работе [50]. В ней предлагается сводить неодносвязную область к односвязной введением дополнительных разрезов «методами теории функций комплексного переменного». С одной стороны, действительно, всегда можно добавить в тело число разрезов, равное количеству внутренних полостей для получения односвязности. С другой стороны, в классической математике и, в том числе, в теории функций комплексного переменного доказывается лишь существование решения, но не исследуется количество операций, необходимых для его построения, и тем более не рассматриваются проблемы конечной точности работы компьютера с действительными числами. Разрез, являющийся отрезком прямой линии, между точками двух несвязных контуров может быть введен только в том случае, если он не пересекает никаких контуров в задаче. Следовательно, для каждой такой пары точек необходимо проверить пересечение будущего разреза с каждым отрезком существующих контуров, в случае обнаружения пересечения взять другую пару точек и так далее. Если обозначить при помощи Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ n суммарное число вершин (равное числу отрезков) в исходных замкнутых контурах триангулируемой области, то проверка пересечения некоторого отрезка со всеми контурами тела потребует O(n) операций, а выбор пары концевых точек разреза возможен O(n2 ) способами. Полное количество необходимых разрезов на единицу меньше числа контуров может быть O(n). Таким образом, непосредственная реализация такого метода устранения неодносвязности будет характеризоваться временем выполнения O(n4 ). Для тысяч входных контуров и десятков тысяч вершин на них не стоит ожидать завершение этого алгоритма в разумное время на современных компьютерах. Но даже если дождаться его завершения, результат может оказаться неудовлетворительным: введенные разрезы могут касаться существовавших контуров. Это происходит от того, что тест пересечения двух отрезков подвержен ошибкам из-за конечной точности оперирования компьютером с действительными числами, и сам по себе является сложной вычислительной задачей [51]. После получения односвязной области в статье [50] предлагается триангулировать ее, образуя треугольники из пары последовательных ребер с контура границы и нового ребра, если оно не пересечет контур в других местах. Здесь уместны ровно те же замечания касательно времени работы и потенциальной опасности ошибок округления. Поэтому в данной работе использовался более сложный в реализации алгоритм триангуляции многосвязных невыпуклых областей, который, однако, обеспечивает лучшую производительность и на практике меньше «страдает» от ошибок округления. В работе [42] предлагается алгоритм, основанный на разбиении области на ряд монотонных многоугольников, которые пересекаются любой прямой параллельной некоторой заданной не более чем в паре точек. Это свойство позволяет триангулировать монотонные многоугольники за линейное время относительно сложности их границы. Алгоритм, приведенный далее, использует родственные идеи и имеет Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ ту же оценку количества операций по порядку величины, однако не выполняет явно разбиение исходной области на монотонные куски, что упрощает реализацию и дает некоторый выигрыш во времени. Алгоритм триангуляции Алгоритм во время своей работы оперирует со структурой данных, введенной в предшествующем разделе. В начале работы алгоритма эта структура уже не пуста и содержит вершины и ребра с контуров области, которую необходимо триангулировать. Каждый контур указывается одним из полу-ребер, входящим в контур и оставляющим внутренность области справа. По исходному полу-ребру легко перечислить все оставшиеся вершины и так же ориентированные полу-ребру контура за линейное время относительно сложности контура. Все индексы полу-ребер сортируются согласно положению вершины, из которой они исходят, в лексикографическом порядке: сначала по x, затем при равных x по y. В алгоритме сортируются полу-ребра вместо вершин, потому что полу-ребро однозначно определяет и вершину, и проходящую через нее границу. Пример, когда через одну вершину проходит две границы приведен на рис. 2.4. Следующим этапом является выборка полу-ребер в отсортированном порядке и совершения действий в зависимости от местного характера границы, заключающихся в построении треугольников в окрестности текущей вершины, что приводит к постоянному изменению границ оставшейся области, а также в поддержании дополнительной структуры данных. Можно представить, что обработка сопровождается движением вертикальной прямой слева направо, при прохождении которой через очередную вершину производится некоторая работа (разновидность plane sweep algorithm). Вместе с движением прямой поддерживается сбалансированное бинарное дерево полу-ребер, которые ей пересекаются и оставляют внутренность об Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ ласти выше (рис. 2.1). Предикатом, отвечающим за размещение элементов в дереве, является сравнение их y-координаты. В начале алгоритма и в его конце дерево является пустым. Количество пересекаемых прямой ребер всегда четное, они распадаются на пары, соответствующие верхней и нижней границе одной из частей области. Поэтому в дереве хранится лишь нижнее полу-ребро пары. С каждым элементом в дереве также ассоциировано хранится ссылка на вершину (определяемая индексом соответствующего полу-ребра), которая была пройдена вертикальной линией в текущей части дерева последней. Эта вершина принадлежит либо одному из пересекаемых ребер, либо является пройденной объединяющей вершиной (см. далее). Следующее утверждение является истинным (инвариант цикла) перед и после выполнения действий над любой вершиной: часть текущей границы, состоящей только из ребер, целиком лежащих левее вертикальной линии, имеет внутренние углы только большие или равные (является выпуклой направо).

a b a end vertex A e0,v0 v3 e3 v0 v2 e2 e0 e3,v b start vertex B e2,v2 e1,v b a split vertex a b merge vertex ba v1 e ab regular vertices a b A B Рис. 2.1. Деревья пересекаемых ребер для двух различных вертикальных линий и классификация типов вершин в алгоритме триангуляции.

Работа, производимая при прохождении вершины зависит от ее типа, который в свою очередь определяется направлением двух граничных полу-ребер a и b, выходящих из нее: a оставляет триангулируемую область Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ справа, b — слева (рис. 2.1). Для упрощения изложения предположим, что a, b не являются строго вертикальными, хотя в программе такого предположения не делалось (подобные случаи встречались и корректно обрабатывались). После выполнения работы a, b добавляются в дерево пересекаемых движущейся прямой полу-ребер, если они оставляют триангулируемую область выше и направлены в правую сторону, либо исключаются, если — в левую. Классификация вершин была заимствована из [42] и имеет непосредственное отношение к разбиению многоугольника на монотонные части. При прохождении начальной (start) или расщепляющей (split) вершины число ребер, пересекаемых вертикальной линией, увеличивается на два. Наоборот, конечная (end) или объединяющая (merge) вершины сокращают это число на два. При прохождении обычной (regular) вершины, число пересекаемых ребер не изменяется. Обработка начальной вершины (ax > 0, bx > 0, [a, b] < 0) наиболее проста. Никаких действий, за исключением запоминания нижнего пересекаемого ребра и самой вершины в дереве при этом не производится.

b

p Рис. 2.2. Обычная вершина границы. Последовательное создание треугольников определяется углом между предшествующей парой ребер границы. Серым обозначена область, подлежащая триангуляции.

В обычной вершине (ax bx < 0) внимание обращается на предшествующую последовательную пару ребер границы (рис. 2.2): если они оба пройдены вертикальной линией и образуют между собой угол меньший, то они становятся ребрами нового треугольника, а третье ребро треугольника заменяет их на границе, после чего проверка возможности создания еще одного треугольника повторяется снова. Текущая вершина запомина Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ ется в дереве пересекаемых ребер, как последняя обработанная для данного участка.

b a b a Рис. 2.3. Триангуляция в окрестности объединяющей (слева) и конечной (справа) вершин сводится к процедуре, разработанной для обычных вершин.

Конечную вершину (ax < 0, bx < 0, [a, b] < 0) можно не выделять от обычных, разве что не выполнять проверку на угол вводимых треугольников из-за того, что она всегда будет успешно пройдена. Обработка объединяющей вершины (ax < 0, bx < 0, [a, b] > 0) эквивалентна действием над парой обычных вершин для участка границы выше и ниже (рис. 2.3). Из двух пересекаемых вертикальной линией донных ребер в окрестности вершины одно удаляется из дерева, а текущая вершина ассоциируется с оставшимся в дереве ребром.

v b a bottom Рис. 2.4. Обработка расщепляющей вершины: слева — ситуация до обработки, показаны порождаемые треугольники, справа — ситуация после. Серым цветом выделена та область тела, которая еще не триангулирована. Bottom — ближайшее снизу полу-ребро, пересекаемое вертикальной линией, v — ассоциированная с ней вершина.

Действия, производимые в расщепляющей вершине (ax > 0, bx > 0, [a, b] > 0), показаны на рис. 2.4. В дереве находится ближайшее к вершине Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ нижнее (bottom) полу-ребро, с которым хранится ссылка на последнюю обработанную на этом участке вершину v. Эта вершина «видна» из любой точки, лежащей на отрезке вертикальной линии между нижней и верхней границей текущей части многоугольника. Поэтому расщепляющая вершина соединяется новым ребром с v. Также кроме v из расщепляющей вершины могут быть видны соседние (и только соседние в силу выпуклости) с ней ребра текущей границы. Эти ребра и расщепляющая вершина служат основой для вводимых треугольников. Операция сортировки имеет трудоемкость O(n log n). При проходе каждой вершины в сбалансированном бинарном дереве проводится поиск нижнего полу-ребра, после чего, возможно, производится обновление ассоциированных данных ребра, его удаление, замена или добавление ребра, лежащего непосредственно выше. Это требует O(log n) операций для обработки одной вершины [52], или O(n log n) — для всех вершин. Количество прочих действий, совершаемых в каждой вершине пропорционально количеству создаваемых треугольников. А так как общее количество треугольников в триангуляции есть O(n), то так же ограничено и полное количество действий, не связанных с деревом пересекаемых ребер. Таким образом, сложность алгоритма есть O(n log n).

2.3.

Оптимальная триангуляция Делоне Триангуляцией Делоне (Delaunay) [42] набора точек называется та кая триангуляция, которая максимизирует в лексикографическом смысле набор величин всех углов треугольников, отсортированный по возрастанию (т. е. в первую очередь максимизируется минимальный угол). Критерием триангуляции Делоне является отсутствие внутри описанной окружности любого треугольника каких-либо других вершин (для ограниченной триангуляции рассматривается только часть окружности, Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ лежащая по ту же сторону границы, что и треугольник). Если вершины отсутствуют и на самой окружности, то триангуляция Делоне заданного набора точек определяется однозначно.

Рис. 2.5. Два принципиальных сценария отношения описанных окружностей и вершин смежных непересекающихся треугольников.

Рассмотрим пару смежных непересекающихся треугольников (пересечение совпадает с общим ребром треугольников). Легко видеть, что возможны лишь два сценария: либо обе описанные окружности треугольников не содержат всех вершин другого треугольника, либо обе они содержат все четыре вершины (рис. 2.5).

a a b b a+b>p a+b

На практике триангуляция Делоне обычно строится из произвольной триангуляции при помощи ряда операций «перещелкивания» (ip) ребра между парой смежных треугольников, для которых не выполнен критерий (рис. 2.6). До перещелкивания обе описанные окружности включают все вершины смежных треугольников. После перещелкивания описанные вокруг новых треугольников окружности по определению имеют минимальный радиус из всех окружностей, включающих все три вершины треуголь Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ ника. Следовательно, после перещелкивания радиусы новых описанных окружностей становятся меньше обоих радиусов окружностей до перещелкивания. По теореме синусов синус угла равен отношению противолежащей стороны к удвоенному радиусу описанной окружности. Поэтому углы треугольников, расположенные против внешних ребер условного четырехугольника, (один из которых имел минимальную величину до изменения триангуляции) в результате перещелкивания увеличиваются. Углы против общего ребра треугольников (диагонали) равны сумме двух углов в другой конфигурации. Следовательно, операция перещелкивания для треугольников с невыполненным критерием Делоне увеличивает минимальный угол. Алгоритм построения оптимальной триангуляции Делоне сводится к проверке всех внутренних ребер триангуляции и смежных с ними треугольников, а также ребер треугольников, порожденных в результате перещелкивания. Все вершины двух смежных треугольников лежат на одной окружности, если сумма противолежащих общему ребру (хорде окружности) углов равна. Если эта сумма больше, то четвертая вершина попадает во внутренность окружности. Таким образом, критерием является выражение (рис. 2.6) +. Непосредственное вычисление углов из координат вершин требует обращения к медленным обратным тригонометрическим функциям. Чтобы этого избежать, рассмотрим эквивалентное выражение sin cos + sin cos = sin( + ) 0. Синусы и косинусы углов равны отношению векторных и скалярных произведений векторов ребер треугольника, отнесенные к произведению длин этих векторов. Но так как сравнение ведется с нулем, то для проверки достаточно вычисления лишь простых скалярных и векторных произведений. Дополнительного выигрыша в скорости можно достичь, заметив, что Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ синусы углов треугольника всегда положительны. Поэтому если cos 0 и cos 0, то критерий выполнен и, наоборот, если cos < 0 и cos < 0, то критерий нарушен. Проверка критерия существенно затрудняется конечной точностью вычислений в компьютере и наличием ошибок округления. Например, когда sin cos и sin cos примерно равны по модулю и противоположны по знаку, то знак неравенства может быть определен ошибочно. Это может привести к тому, что алгоритм зациклится и будет выполнять перещелкивания в соседних парах треугольников до бесконечности. Чтобы избежать этого, вычисления можно проводить с большей точностью, что отрицательно скажется на производительности. В статье [51] исследуется этот и другие геометрические предикаты и показывается, как избежать ошибок и существенных потерь производительности, используя так называемую адаптивную точность вычислений. В данной работе для избежания зацикливания использовалась проверка sin cos + sin cos >. ( > 0) В этом случае критерий Делоне может нарушаться, но только для пар треугольников, описанные окружности которых практически совпадают, что приводит к также достаточно качественной триангуляции.

2.4.

Поддержание заданной плотности сетки Автором был разработан алгоритм построения равномерной треуголь ной сетки, который исходит из заданной триангуляции тела, внося в нее минимальные изменения. В качестве исходной обычно берется триангуляция тела с предшествующего шага интегрирования, вершины которой смещены на произведение скорости тела в этих точках и величины шага интегрирования. Можно так определить порог на временной интервал между шагами (см. раздел 2.9), что даже после смещения триангуляция остается Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ невырожденной (никакие два треугольника не имеют пересечение большее их общего ребра). Мелкость сетки и допустимая степень ее неравномерности опредеi ляются двумя числовыми параметрами алгоритма: lmin > 0 и > 1 со ответственно. Как будет доказано ниже алгоритм гарантирует, что (при отсутствии слишком острых углов с внутренней стороны любой границы области [53]) для всех треугольников в сетке выполняются условия на их минимальный и максимальный размеры (теорема 1). Минимальность вносимых изменений в сетку (появление и удаление вершин) обеспечивается за счет увеличения параметра (в расчетах использовалось = 1.1). Идея заключается в том, что каждое ребро в триангуляции может иметь длину от некоторого минимального размера lmin до размера большего в 2 раз. При этом допустимые пределы на размеры объектов сетки увеличиваются за счет ее равномерности. В результате растяжения тела вводятся новые ребра длиной не меньше lmin. Если бы размер максимально допустимого ребра не превосходил более чем вдвое размер минимального ребра, то введение дополнительных вершин могло порождать слишком короткие ребра. Когда растяжение тела сменится сжатием (например, в случае малых колебаний тела), то у нового ребра будет некоторый запас уменьшения длины, прежде чем его потребуется удалить из триангуляции. Отличие предлагаемого алгоритма от уже известных способов построения сетки [53] заключается в ориентации на расчет эволюционных процессов, что подразумевает создание сетки близкой к равномерной. Также при расчете по явным сеточно-характеристическим схемам оказывается важным наличие нижней границы на высоту в любом треугольнике сетки (см. главу 5). Алгоритм поддержания заданной плотности сетки состоит из следующих шагов:

Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ 1. Построение триангуляции Делоне из данной изначально и поддержание этого свойства в дальнейшем. 2. Сокращение граничных вершин до тех пор, пока имеются граничные b i ребра короче lmin < lmin.

3. Введение дополнительных вершин в центрах граничных ребер длиннее b b lmax = 2lmin. i 4. Сокращение внутренних вершин, одно из ребер которых короче lmin.

5. Введение новых внутренних вершин в центрах описанных окружноi i стей треугольников сетки, если их радиус больше Rmax = lmin, а центр попадает внутрь тела. 6. Проверка границы триангуляции на наличие в ней треугольника, сразу два ребра которого являются граничными, а высота меньше hb. В min случае обнаружения таких треугольников в их окрестности вводятся (и, возможно, удаляются) дополнительные внутренние вершины, чтобы обновленная триангуляции Делоне более не включала подобных вырожденных треугольников. Практика показывает, что вероятность их возникновения крайне низка и данный шаг алгоритма можно при желании опускать.

b b i Величины lmin, lmax, hb являются функциями lmin, и будут опреmin делены ниже. Рассмотрим наиболее сложные не описанные выше шаги алгоритма подробнее. 2.4.1. Сокращение граничных вершин Простейшая реализация шага 2 алгоритма не представляет никаких затруднений. Необходимо лишь удалять одну из вершин слишком коротких граничных ребер до тех пор, пока на границе тела не останется ни одного Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ короткого ребра. Недостаток такого пути проявляется в расчете контактного взаимодействия нескольких тел (глава 3). Дело в том, что в результате удаления вершин на вогнутой границе площадь тела может немного увеличится, и оно «залезет» на контактирующий с ним объект.

?

Рис. 2.7. Необходимые действия для сокращения граничного ребра без увеличения площади тела в зависимости от выпуклости или вогнутости тела в окрестности. Черными точками показаны существующие вершины, белыми точками — добавляемые, крестиками — удаление вершин или отдельных ребер.

Чтобы избежать увеличения площади тела в следствие сокращения граничных вершин предлагается исследовать окрестность короткого граничного ребра и действовать одним из перечисленных способов. Если в районе одной из граничных вершин тело выпукло, а в районе другой — вогнуто (рис. 2.7, а), то удаляется «выпуклая» вершина. Если тело выпукло в окрестности обоих вершин (рис. 2.7, б), то они удаляются, но, чтобы сокращение площади не было чрезмерным, вводится новая вершина в центре бывшего граничного ребра. Когда тело вогнуто, то для принятия решения необходимо определить, попадает ли точка пересечения продленных двух соседних граничных ребер внутрь треугольника, гранью которого является короткое граничное ребро. Если да (рис. 2.7, в), то также удаляются обе вершины и вводится одна в исследованной точке пересечения. Если нет (рис. 2.7, г), то удаляется лишь только короткое ребро и содержащий его треугольник, а два прочих ребра удаленного треугольника становятся граничными. С другой стороны, сокращение площади при выравнивании плотности сетки на границе не компенсируется другими процессами, и при моделировании продолжительных динамических задач тела могут стать заметно меньше в размерах. Если такая трудность возникает, то необходимо Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ уменьшить мелкость сетки, и тогда потери площади также сократятся.

2.5.

Обоснование корректности алгоритма Разобьем анализ алгоритма на формулировку и доказательство ряда утверждений. Утверждение 1. В триангуляции Делоне вообще и в построенной сетке, в частности, у любой вершины всегда есть ребро, соединяющее ее с наиболее близкой другой вершиной.

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

Доказательство. Все точки, расположенные к данной вершине ближе, чем длина ее самого короткого ребра, находятся внутри объединения описанных окружностей вокруг соседних с вершиной треугольников (рис. 2.8). Применяя критерий триангуляции Делоне, приходим к сформулированному утверждению. Утверждение 2. После выполнения шага 4 алгоритма расстояние между любой вершиной сетки, не лежащей на границе, и любой другой вершиi ной сетки не менее lmin. Это свойство поддерживается на последующих шагах.

Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ Доказательство. Шаг 4 алгоритма устраняет в сетке все внутренние верi шины, имеющие ребра длиннее lmin. Следовательно, длина самого коротi кого ребра у внутренних вершин не меньше lmin. Доказательство вытекает из утверждения 1. Это утверждение можно было бы перенести на длину любого внутреннего ребра, если бы не существовало внутренних ребер, обе вершины которых принадлежат границе. Справедливость утверждения в таком случае зависит от состоятельности контуров тел, поданных на вход алгоритма. Необходимо, чтобы контуры не приближались друг к другу на расстояние i менее lmin, а также не имели слишком острых внутренних углов (подробнее в подразделе 2.7.1). При тех же условиях расстояние от любой граничной b вершины до ближайшего своего соседа не меньше lmin.

Далее будем называть внутренним треугольником сетки такой треугольник, не более одной вершины которого принадлежит границе тела.

b i Утверждение 3. Если lmax 2lmin, то при выполнении шага 5 алгоритма центр описанной окружности любого внутреннего треугольника лежит внутри тела.

limin j2 j3 j4 j1 R lbmin l min i lbmin j j2 j3 j R Рис. 2.9. Внутренний треугольник сетки, центр описанной окружности которого лежит за пределами тела (обозначено серым цветом). Никакие вершины сетки не могут располагаться внутри представленных окружностей по свойствам триангуляции Делоне. Размер и положение минимального граничного ребра зависит от i :

1 2 i i > (слева) и 1 i i (справа).

Доказательство. Предположим противное: в текущем состоянии сетки, являющейся триангуляцией Делоне, нашелся внутренний треугольник, центр Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ описанной окружности которого лежит за пределами тела (рис. 2.9). По свойствам триангуляции Делоне никакие другие вершины сетки не могут находиться внутри этой описанной окружности, а также внутри окружноi b стей с центрами в вершинах треугольника и радиусами lmin или lmin в зави симости от типа вершины. Это накладывает ограничение снизу на размер граничного ребра, которое пересекается описанной окружностью треугольника. Если радиус описанной окружности достаточно большой, то имеет место ситуация, представленная на рис. 2.9 справа, и наиболее короткое граничное ребро имеет концы в точках пересечения указанных окружностей. Углы секторов, образованных линиями, соединяющими вершины треугольника и эти точки с центром описанной окружности, обозначены как i. Критерием «достаточно большого радиуса» является условие 1 2 i i i. Если учесть, что расстояние между точками не меньше lmin 2 b для внутренних узлов и lmin для граничных, то i lmin 1,2,3, sin 2 2R b 4 lmin sin. 2 2R Условие на длину пересекаемого граничного ребра можно записать следующим образом l 2R · b sin( 1 1, i i ), если если 1 2 1 i i i i >, 2. Минимум выражения lb = 2R достигается, когда треугольник имеет такой радиус, что i lmin 1,2,3 =, sin 2 2R b 4 lmin sin =, 2 2R 1 i = i. Следовательно, i b lmin lmin 3 arcsin + arcsin = 2R 2R i lmin < sin 2R b i lmax lb = 2R > 2lmin.

Получено противоречие с условиями утверждения.

Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ Достаточное условие, фигурирующее в условии утверждения, может быть ослаблено, но для данного исследования это не потребуется. Утверждение 4. Выполнение шага 5 алгоритма всегда завершается, и b i если lmax 2lmin, то длины ребер всех внутренних треугольников после i завершения не превышают 2Rmax.

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

b i Если lmax 2lmin, то согласно утверждению 3 не может найтись внутреннего треугольника, описанная окружность которого имеет центр за пределами тела. Поэтому после завершения шага 5 алгоритма не остается внутренних треугольников, радиус описанной окружности которых i более Rmax, и, следовательно, ни одно ребро внутренних треугольников не i длиннее 2Rmax.

Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ 2.6.

Размеры внутренних треугольников сетки b i Теорема 1. После окончания работы алгоритма при условии lmax 2lmin справедливы следующие оценки для внутренних треугольников сетки: 1i i lmin Ri lmin, i lmin 1 arcsin 6 2 2 3= i max 21 4 i (lmin )2 2 4 i 3lmin i li 2lmin, i 2 arcsin ni 2 1, 2 3 (2.1) 2 12, i min 332i 2 S i (lmin ), 4 i pi 3 3lmin, + 1 i 42 1 lmin, 1i lmin hi 42 1 i i lmin ri lmin. 2 42 ( 42 1 + 23 ) где R — радиус описанной окружности вокруг треугольника, l — длина ребра, — угол треугольника, n — число смежных треугольников у вершины, S — площадь треугольника, p — периметр треугольника, h — высота, опущенная из произвольной вершины треугольника, r — радиус вписанной окружности. Причем эти оценки являются точными и могут быть достигнуты. Доказательство. Оценки для длины ребра и ограничение на радиус описанной окружности сверху были получены в утверждениях 2, 4. Минимальный радиус описанной окружности для внутреннего треугольника сетки, очевидно, достигается для правильного треугольника, i все ребра которого имеют минимально возможную длину lmin. Записывая теорему косинусов для «трети» этого правильного треугольника, получаi i ем уравнение на радиус: (lmin )2 = 2(Rmin )2 (1 cos 2 ), откуда находим 3 i Rmin = 1 li. 3 min Следует отметить, что треугольник, вписанный в окруж Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ ность минимального радиуса, — единственен (с точностью до произвольного сдвига и поворота), тогда как с увеличением радиуса количество форм вписанных треугольников возрастает. Поэтому экстремальные характеристики треугольников в большинстве случаев следует искать либо у единственного минимального правильного треугольника, либо у одного из треугольников, вписанных в окружность максимального радиуса.

l min i j min i j max i i lm j min i in R ma i x Рис. 2.10. Наиболее вырожденный внутренний треугольник сетки: два его угла равны минимальному, третий — максимальному. Этот треугольник также имеет минимальную площадь, высоту и радиус вписанной окружности.

Возможные углы в треугольнике сетки можно оценить из теоремы синусов, согласно которой отношение стороны треугольника к синусу противолежащего угла равно диаметру описанной окружности. Следовательно, минимальный угол достигается в треугольнике, вписанном в окружi ность максимального радиуса Rmax, против стороны минимальной длины i lmin : sin i = min i lmin i 2Rmax = 1 2.

На рис. 2.10 приведен треугольник, имеющий сразу два минимальных угла. В нем же достигается и максимальный угол 1 i max = 2 arcsin 2. Для использованного в данной работе значения = 1.1 получается i 27, i 126. max min Минимальные и максимальные углы треугольников непосредственно отвечают на вопрос возможного числа смежных треугольников у одной вершины: ni = min условиям Делоне. Определение разброса площадей треугольников сетки. Площадь треугольника есть половина произведения длины одной его стороны a на высоту h, опущенную из противолежащей стороне вершины. Среди треуголь2 i max, ni = max 2 i min. Рис. 2.11 демонстрирует, что кон фигурации с таким количеством смежных треугольников удовлетворяют Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ i lm mi ax in j min i jimax R ma i x Рис. 2.11. Минимальное и максимальное число смежных треугольников у одной вершины. Описанные окружности не включают вершины других треугольников, поэтому данное расположение треугольников допустимо в триангуляции Делоне.

R h a Рис. 2.12. Для вписанных в одну окружность треугольников максимальная площадь будут у одного из равнобедренных треугольников.

ников, вписанных в окружность радиуса R, максимальную площадь имеет один из равнобедренных треугольников (рис. 2.12): a = 2 R2 (h R)2 = 2 2hR h2, 1 S = ah = 2h3 R h4, 2 S 2 = 6h2 R 4h3. h Приравнивая производную от площади к нулю, узнаем, что максимальная площадь достигается при h = 3 R, a = 3R, т. е. в правильном равносто2 роннем треугольнике, i Smax = 33 i 2 4 (Rmax ) = 332i 2 4 (lmin ).

Минимальная площадь может быть либо у единственного треугольi ника, вписанного в окружность минимального радиуса Rmin : S1 = 3i 2 4 (lmin ), R 33 i 2 4 (Rmin ) = либо у наиболее вырожденного треугольника, вписанного в 1 2 i sin i (lmin )2. Сиmax окружность максимального радиуса (рис. 2.10): S2 = Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ нус максимального угла можно выписать явно: sin i = sin( 2 arcsin max 1 1 ) = sin(2 arcsin ) = 2 2 1 1 1 1 1 2. = 2 sin(arcsin ) cos(arcsin ) = 2 2 При = 1 обе площади совпадают: S1 = S2, при > 1 — меньше S2. Минимальный периметр треугольника, очевидно, равен трем миниi мальным длинам ребра: pi = 3lmin. Для определения максимального пеmin риметра рассмотрим треугольник с углами,,. Согласно теореме синусов можем записать (, ) = arg max p = arg max(,, 1 1 1 + + ). sin sin sin( ) Для отыскания максимума приравняем производные выражения по и нулю: cos( ) cos = 0, sin2 sin2 ( ) cos( ) cos = 0. 2 sin sin2 ( ) Функция cos sin = cos убывает на интервале (0;

), поэтому sin2 cos следует равенство аргументов = = sin из равенства функций =. Таким образом, максимальный периметр достигается в правильном треугольнике, i i вписанным в максимальную окружность: pi = 3 3Rmax = 3 3lmin. max Высота треугольника равна произведению длины одной из его сторон на синус прилежащего к ней угла, который можно выразить, используя теорему синусов, через длину противоположной стороны и радиус описанной окружности: l1 l2. 2R Поэтому минимальная высота будет в треугольнике, имеющем два наиh= меньших ребра и вписанном в окружность максимального радиуса (рис. 2.10): hi = min 1i 2 lmin.

Наибольшая высота достигается в равнобедренном треуголь нике, имеющем форму иглы, в котором основание имеет минимальную Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ длину. Используем формулы, полученные при исследовании площади треугольников:

i i lmin = 2 2hi Rmax (hi )2, max max i hi = Rmax + max 1 i i 4(Rmax )2 (lmin )2 = + 1 i 42 1 lmin.

Радиус вписанной в треугольник окружности равен отношению удвоенной его площади к периметру. Очевидно, что максимальным радиус будет у правильного треугольника, вписанного в наибольшую окружность: 2S i i i rmax = max = lmin. i 2 3 3Rmax Минимальный радиус вписанной окружности достигается в треугольнике с рис. 2.10:

i rmin = i 2Smin i 2Rmax i 2lmin + sin i max = ( 42 1 i 2 22 (lmin ) 43 i + 42 1 )lmin 42 1 = li. 2 ( 42 1 + 23 ) min Некоторые оценки этой теоремы были также экспериментально проверены вычислением минимальных и максимальных характеристик треугольников порожденных сеток.

2.7.

Допустимые размеры длины граничного ребра Рассмотрим в данном разделе такие треугольники, одно из ребер ко торых является граничным. Размеры (углы и т. д.) этих треугольников будем снабжать верхним индексом b. Если на границе области построения сетки разрешить тот же диапаi зон длин для ребер треугольников, т. е. допустить lb > 2lmin, то возникает ситуация, когда внутренний узел сетки может приблизиться к ее границе сколь угодно близко (рис. 2.13, слева). Минимальное расстояние от внутреннего узла до границы обуславливает шаг интегрирования (см. главу 5) Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ при решении динамических уравнений, поэтому сокращение расстояния, тем более неконтролируемое сокращение, крайне нежелательно. Это с неизбежностью приводит к тому, что пределы размеров для граничных ребер должны быть снижены так, чтобы внутренний узел всегда находился на известном минимальном удалении от границы. Отношение максимального и минимального размеров граничных ребер необходимо оставить таким же (2), как для внутренних ребер по причинам, изложенным в разделе 2.4.

l min i h lbmax l i min >R max i l min i b j max l min b i j min j min h i R max b Рис. 2.13. Две причины уменьшения высоты приграничного треугольника по сравнению с минимальной высотой внутренних треугольников hi : превышение радиусом описанной окружmin i ности, центр которой лежит вне тела, порога Rmax (слева) и уменьшение длины граничного ребра (справа).

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

b b В данной работе было решено выбрать пределы lmin, lmax таким об разом, чтобы минимальные высоты в обоих случаях, изображенных на рис. 2.13 были равны. Из ситуации, изображенной слева: 1b i (hb )2 = (lmin )2 (lmax )2. min 4 Из ситуации, изображенной справа: hb min i b b b lmax lmin lmin lmin = = =. i 2Rmax 2 Подставляя выражения одно в другое, получаем hb min 1 i = lmin, 1 + b lmin 2 i = lmin, 1 + b lmax 42 i = lmin. 1 + b Заметим, что величина lmax при любом удовлетворяет условию теоремы 1.

Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ Следовательно, если от любой вершины в сетке сместиться на расстояние не большее hb (только внутрь области для граничных вершин), min то конечная точка окажется внутри одного из треугольников, смежных с исходной вершиной. Для треугольников с двумя граничными ребрами это непосредственно проверяется последним шагом алгоритма. 2.7.1. Минимальный угол границы тела l i lmin j l b min b Рис. 2.14. Минимальный угол границы тела, при котором выполняется условие на допустимый размер внутреннего ребра.

Любой алгоритм построения треугольной сетки налагает определенное ограничение на минимальный угол области, в которой строится сетка [53]. Не является исключением и предлагаемый в этом разделе алгоритм. Если рассмотреть треугольник со сторонами a, b, c и углом против стороны a, то справедливо равенство (теорема косинусов): a2 = b2 + c2 2bc cos. Применяя это равенство к случаю, изображенному на рис. 2.14, приходим к неравенству 1 + 44 b i b lmin = lmin li = 2lmin (1 cos ) 2 необходимо > 66. cos 1 + 44. Поскольку > 1, то > 63.8. Для использовавшегося в расчетах = 1. Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ 2.7.2.

Обработка треугольников с двумя граничными ребрами Выше было показано, что для внутренних треугольников сетки и треугольников, одно из ребер которых принадлежит границе, можно гарантировать, что любая их высота будет не меньше величины hb. Однако для min треугольников, имеющих сразу два граничных ребра этот порог должен быть понижен. Однако есть возможность не изменять оценку высоты, а проверять все подобные треугольники, и в случае обнаружения нарушения оценки, вносить соответствующие изменения в сетку. Именно это имеем место на 6-м шаге алгоритма.

? Rimax Rimax h

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

На рис. 2.15 слева представлен треугольник, для которого нарушается оценка высоты. Если удастся ввести внутри его описанной окружности новую внутреннюю вершину так, что она будет удалена не менее чем i на Rmax от вершин треугольника, то в обновленной триангуляции Делоне проблемный треугольник исчезнет (рис. 2.15 справа). При этом потребуi ется удалить все внутренние вершины, расположенные ближе чем lmin от вновь вводимой, а также повторить шаг 5 алгоритма в окрестности, чтобы восстановить возможно нарушенную введением вершины сетку. Рассмотрим, как часто возникают подобные вырожденные треуголь Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ ники и критерии возможности исправления сетки. Оговоримся, что будем исследовать вопрос только для разумных очертаний тела, в частности, без резких изломов (см. подраздел 2.7.1) и областей утончения до величин поi рядка lmin.

Утверждение 5. Если треугольник с двумя граничными ребрами имеет слишком маленькую высоту: h < hb, то оба его граничных ребра короmin i че lmin, а приграничный угол треугольника (угол излома границы в этом месте) лежит в диапазоне b hb 1 lmin 1 min = 2 arccos b 2 arcsin i = 2 arcsin 2 arccos. 2 2Rmax lmin 1 + Доказательство. Поскольку мы считаем, что очертания тела лишены чересчур острых углов, то исследуемый треугольник должен иметь центр описанной окружности внутри тела, и, следовательно ее радиус не превосi ходит Rmax. В начале раздела было показано, что в любом треугольнике, вписанном в такую или меньшую окружность и имеющем лишь одну стоi рону короче lmin, высота всегда удовлетворяет условию h hb. Следоваmin i тельно, чтобы h < hb, оба граничных ребра должны быть короче lmin. min Наибольший угол достигается, когда треугольник вписан в окружi ность максимального радиуса Rmax, и в нем имеется два минимальных угла b против сторон длиной lmin :

max b lmin = 2 arcsin i. 2Rmax Наименьший угол треугольника, когда еще h hb, имеет место min b также при длине двух граничных ребер lmin :

min hb cos = bmin. 2 lmin Если считать размеры всех граничных ребер независимыми равноb b мерно распределенными случайными величинами на отрезке [lmin ;

lmax ], то Глава 2. ПОСТРОЕНИЕ НЕРЕГУЛЯРНОЙ ТРЕУГОЛЬНОЙ СЕТКИ вероятность выполнения необходимого условия на ребра вырожденного треугольника для = 1.1 менее 3%. К тому же, следует еще произвести умножение на вероятность угла границы также быть в диапазоне [126 ;

Pages:     || 2 | 3 |



© 2011 www.dissers.ru - «Бесплатная электронная библиотека»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.