Численное моделирование деформационных динамических процессов в средах со сложной структурой

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

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

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

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

Краткие обзоры-классификации используемых численных методов и видов сеток можно найти в работах [5,6].

0.1. Способы дискретизации уравнений механики.

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

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

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

Метод конечных разностей (finite-difference method) основывается на. замене производных в поставленных дифференциальных уравнениях на конечные, недифференциальные комбинации (в простейшем случае — разности) значений, определенных в точках сетки. Приближенное решение состоит из значений U^ (в двумерном случае), представляющих поточечное приближение к решению в узлах сетки (Xi, yj) во время tn:

Ujj ~ и (х{, yjt tn).

Метод конечных объемов (finite-volume method), с другой стороны, основывается на приближении средних значений и в ячейках сетки.

ГУ] +½ r*i+l/2 u (x, y, tn)dxdy.

Jyj-l/2 Jxi-l/2.

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

Программы, реализующие как метод конечных разностей, так и метод конечных объемов, могут хранить рассчитанные параметры Ufj в виде прямоугольной таблицы. В первом случае в таблицу попадут значения в О некоторых точках тела, во втором — усредненные значения по ячейкам.

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

Но алгоритмы численных методов, записанные как действия над элементами таблицы, в обоих случаях будут очень похожи (рис. 1) и могут даже совпадать для линейных уравнений.

Метод конечных элементов (finite element method), примененный, например, в [7−9], относится к вариационным методам [10]. Решение в нем ищется среди конечного множества функций, определенных согласованно с текущей сеткой. Метод конечных элементов наиболее часто используется в случае нерегулярности сетки из-за простоты его обобщения на этот случай.

Методы конечных объемов и конечных элементов высоких порядков аппроксимации с упрощенной, аналитической реконструкцией для неструктурированных сеток получили названия методов спектральных объемов (spectral volume method) [11] и спектральных элементов (spectral element method) [12] соответственно. Некоторые соображения, лежащие в их основе, были заимствованы в данной работе для построения сеточно-характе-ристического метода второго порядка на неструктурированных сетках.

Методы частиц [13] и гладких частиц (smoothed particle hydrodynamics) [14−16], или метод свободных точек [17] применяются при моделировании значительных деформаций тел и разрушений в них, приводящих к отколу множества обломков.

Сеточно-характеристический метод (grid-characteristic method) [1834], используемый в данной работе, ближе всего к конечно-разностным методам в том смысле, что в точках сетки хранятся приближенные значения численного решения в этих точках, и в качестве основы используется дифференциальная форма записи уравнений механики. Однако вместо производных, входящих в исходные уравнения, аппроксимации подвергаются производные вдоль характеристических направлений, что обеспечивает большую точность. Также отличительной особенностью сеточно-характеристического метода является удобство вычисления на границе области интегрирования с привлечением ровно того же числа граничных условий, которые необходимы для корректного замыкания системы дифференциальных уравнений.

0.2. Способы построения сетки в области интегриро.

Методы введения сеток для расчета деформируемых тел классифицируются по трем основным признакам:

• наличию связи между положением сеточных линий и границ исследуемых тел,.

• способностью меняться по мере деформирования моделируемых тел и.

• наличию структуры (регулярности) в сетке.

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

0.2.1. Квадратная регулярная сетка вания.

Рис. 2. Квадратная регулярная сетка в окрестности исследуемого объекта.

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

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

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

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

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

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

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

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

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

• малая кривизна сеточных линий,.

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

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

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

0.2.3. Нерегулярная треугольная сетка.

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

Чаще всего используются нерегулярные сетки (рис. 6) с простей- ' шей формой ячейки — треугольником, хотя иногда применяют четырехугольные и смешанные (треугольно-четырехугольные) нерегулярные сетки. Ячейки более сложной формы используются намного реже.

0.2.4. Изменение сетки при деформировании тел.

Для решения эволюционных задач выбираются несколько промежуточных моментов времени, называемых временными слоями, между исходным состоянием и полным временем моделирования. Решение последовательно отыскивается на каждом слое, причем расчетная сетка может также меняться. В совокупности сетки из всех слоев образуют единую сетку в пространстве координат-времени. Численные методы также классифицируются по виду движения узлов между отдельными сетками при переходе к следующему слою.

В механике сплошной среды принято записывать динамические уравнения в одном из двух видов переменных [5]: эйлеровых — фиксированных в пространстве и лагранжевых — движущихся вместе с точками тела. Аналогично выделяют и два основных типа сеток, но кроме них используют и другие [19]. Рассмотрим далее каждый из них отдельно.

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

• Точки лагранжевой сетки смещаются вместе с точками тела, поэтому границы тела всегда совпадают с сеточными линиями. Однако при наличии сдвиговых деформаций в теле, ячейки лагранжевой сетки могут постепенно вырождаться и пересекать друг друга. Численные методы при вырождении сетки обнаруживают неустойчивость и счет приходится прекращать. Поэтому неизбежен подход, называемый лагран-жева сетка с перестройкой, в котором, как только детектируется приближение сетки к вырожденному состоянию, производится построение новой сетки, а значения в новых узлах интерполируются из прежних значений. Процедура интерполяции сама по себе приводит к потере точности решения, поэтому желательно перестраивать сетку как можно реже. В противном случае имеет смысл обратиться к подвижной сетке. Задачи о соударении в лагранжевых переменных рассматривались в работах [37,38].

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

• По-координатная комбинация разных типов движения сетки рассматривается в статье [19]. Например, если правая и левая границы области интегрирования неподвижны, то дг-координата узлов может быть эйлеровой и не храниться в памяти программы, а у-координата — лагранжевой.

• Совместная эйлерово-лагранжева сетка представлена в работе [39]. Подобная сетка имеет неподвижные точки внутри тела и движущиеся вместе с поверхностью тела точки на ее границе. Основная сложность работы с такими сетками — обеспечить сопряжение двух типов узлов в приграничной полосе. В другой работе [40] исследование задач о соударении в эйлеровых переменных приводило к необходимости использовать метода маркеров.

• Иной вариант совмещения эйлерового и лагранжевого подхода предлагается в методе «частиц в ячейке» Ф. Харлоу [13,17]. Метод использует две сетки частиц: лагранжеву сетку для описания частиц, представляющих элементы жидкости, движущихся по неподвижной эйлеровой сетке. Другой пример совмещения неподвижной сетки с подвижными частицами дает работа [41].

• Неструктурированная лагранжева сетка с перестройкой (изменением «соседства») используется в методике «Медуза», также изложенной в [17]. Область интегрирования разбивается на построенные вокруг заданного набора точек многоугольные ячейки («глобулы»), во многом схожие с ячейками диаграммы Вороного [42]. Заметную трудность о в этой методике представляет расчет граничных ячеек, поскольку их стороны не совпадают с контуром области интегрирования.

• В данной работе предлагается численный метод, основанный на нерегулярной лагранжевой сетке с локальной перестройкой. Для достиже-' ния порядков аппроксимации выше первого в этом методе рассчитываются узлы, находящиеся не только в вершинах треугольников сетки, но также и дополнительные узлы, например, лежащие в центре ребер. Поэтому первые узлы являются лагранжевыми, а вторые — подвижными.

Основные результаты и выводы диссертации.

Аналитическим образом произведено спектральное исследование матриц коэффициентов системы уравнений теории упругости, выписанной в произвольной прямолинейной системе координат. В компактной форме получены выражения для всех собственных значений этих матриц Л (1.22), их левых собственные векторов Q (1.23) и векторов взаимного к ним базиса fi-1 (1.24). В предшествующих работах такие выражения были известны только для декартовой системы координат [55,75], а в прочих случаях опре-. делялись численно [19,22,72].

В работе предлагается использовать явное представление сеточно-ха-рактеристических схем, основываясь на произведенном спектральном исследовании, поскольку в записи таких схем входят громоздкие выражения f относительно A, Q, Q-1. В полученной упрощенной записи не требуется решения системы линейных уравнений, обращения и даже перемножения матриц. Полученные выражения приведены к виду, инвариантному относительно размерности пространства (справедливы в 2D и 3D), тогда как запись A, ft, ft'1 отлична для двумерного и трехмерного пространств [75]. В результате чего вместе с упрощением программы достижима более высокая скорость ее работы, а также исключаются численные ошибки, связанные с решением возможно обусловленных систем линейных уравнений.

Для граничных узлов помимо явной записи было предложено использовать двухэтапный метод, причем первый этап не зависит от граничных условий, а второй — от порядка аппроксимации. Разделение на этапы чрезвычайно удобно в программной реализации, поскольку позволяет отдельно отлаживать компактные модули, отвечающие только за перенос значений вдоль характеристик с тем или иным порядком точности либо только за корректировку для конкретного граничного условия. Метод требует решения лишь системы из m-линейных уравнений, где га-число выходящих из области характеристик, тогда как классический подход [19,22] требует решения полной системы из n-уравнений, где п-число переменных в задаче.

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

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

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

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

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

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

Заключение

.