Оглавление
Стр.
Введение 5
Глава 1. Метод сглаженных частиц 33
§1. Основные уравнения динамики жидкости ............................... 33
§2. Основы метода сглаженных частиц...................................... 35
2.1. Усреднение функции но Стсклову....................................... 35
2.2. Первая производная функции........................................... 36
2.3. Функция ядра......................................................... 38
2.4. Дискретизация области................................................ 42
§3. Дискретные уравнения динамики сплошной среды......................... 44
3.1. Плотность и уравнение неразрывности.................................. 44
3.2. Уравнения движения Навье-Стокса...................................... 46
3.3. Вторая производная. Оператор Лапласа................................. 47
3.4. Интегрирование по времени............................................ 49
3.5. Радиус сглаживания................................................... 54
§4. Граничные условия.................................................... 56
4.1. Твердая граница...................................................... 56
4.2. Свободная граница.................................................... 60
4.3. Вычисление гидродинамических нагрузок................................ 61
§5. Алгоритм расчета методом сглаженных частиц .......................... 62
Глава 2. Тестирование метода 63
§1. Решение тестовых задач. Сходимость.................................. 63
1.1. Задача о деформации жидкого эллипса....................... 63
1.2. Задача о ламинарном течении жидкости в плоском канале .... 70
1.3. Задача о течении жидкости по наклонной плоскости........... 75
2
1.4. Задача о падении капли в жидкость............................... 78
§2. Алгоритмы стабилизации вычислений .............................. 81
2.1. Алгоритм регуляризации расчетной сетки......................... 81
2.2. Использование потенциала парного взаимодействия для корректировки свободной поверхности......................................... 87
§3. Задачи на вычисление давления и гидродинамических нагрузок . . 91
3.1. Задача о колебаниях жидкости в прямоугольном бассейне .... 91
3.2. Задача о колебаниях жидкости в прямоугольном бассейне при
наличии режимов обрушения ..................................... 93
§4. Задача о разрушении плотины..................................... 96
4.1. Моделирование турбулентных режимов течения..................... 97
4.2. Постановка задачи .............................................100
4.3. Результаты вычислительных экспериментов........................100
Глава 3. Взаимодействие жидкости с погруженным телом 105
§1. Уравнения движения погруженного тела в жидкости.................105
§2. Граничные условия на поверхности твердого тела..................109
§3. Задача о всплытии плоского цилиндра с круглым основанием в
бассейне с жидкостью.............................................111
3.1. Постановка задачи .............................................111
3.2. Результаты расчетов............................................112
§4. Задача о входе плоского цилиндра в бассейн с жидкостью..........121
4.1. Постановка задачи .............................................121
4.2. Вход в жидкость цилиндров с различными формами оснований.
Результаты расчетов............................................122
§5. Задача о несимметричном входе плоского цилиндра с квадратным
основанием в бассейн с жидкостью................................134
5.1. Постановка задачи .............................................134
5.2. Результаты расчетов............................................135
3
§6. Вход в жидкость плоских круговых цилиндров различной массы . 139
6.1. Постановка задачи ................................................140
6.2. Результаты расчетов...............................................140
Заключение 151
Литераз ура 153
Приложение 1. О погрешности интегральных аппроксимаций 170
Приложение 2. Об аппроксимации полиномов 174
Приложение 3. Об определении положения свободной границы 176
4
(
Введение
Современную науку невозможно себе представить без использования мощнейшего аппарата вычислительной математики, который с каждым годом приобретает все большую популярность. Непрерывный рост производительности компьютерной техники, а также ее доступность для широкого круга исследователей, обуславливает постоянное увеличение общей доли вычислительного эксперимента в научных исследованиях по сравнению с натурными или лабораторными испытаниями. Разумеется, вычислительный эксперимент не претендует на исключительную роль в научных исследованиях, однако, он позволяет значительно снизить потребность в проведении экспериментов, что неизбежно влечет за собой сокращение как временных, гак и денежных затрат на научные исследования.
Рост производительности компьютерной техники сказывается не только на экстенсивном использовании вычислительных технологий (уменьшение шага расчетной сетки, ускорение расчетов благодаря использованию вычислительных мощностей множества ядер или процессоров - параллельное программирование и др.), но и на повышении их качественных характеристик: разрабатываются новые численные методы, постоянным модификациям подвергаются уже известные и широко используемые, расширяется сфера их применения и т.д.
Применение новейших методов позволяет моделировать сложные задачи гидродинамики, включая задачи с большими деформациями свободных и контактных границ между жидкостями и погруженными телами, задачи с межфазными переходами, турбулентные течения, кавитационные обтекания тел, быстропротекающие процессы, такие как подводные взрывы и т.п.
К одному из наиболее сложных для моделирования и представляющих большой практический интерес классов задач гидродинамики относятся задачи со свободными 1раницами. В качестве примеров, где встречаются такого
5
рода задачи, можно привести такие сложные физические процессы как накат волны на наклонный берег, взаимодействие волн с береговыми и донными сооружениями, погруженными в жидкость телами: задачи глиссирования, посадки гидросамолетов на поверхность водоемов и т.д. Практический интерес в таких задачах представляют как кинематические характеристики течений -поле скоростей, положение свободной границы, так и динамические - поле давления, величины гидродинамических нагрузок на твердые границы области течения или погруженные тела.
Исследованию процессов взаимодействия твердых и упругих тел с жидкостью посвящены работы многих отечественных и зарубежных ученых: Л.И. Седова [1], Г.В. Логвиновича [2], В.В. Пухначева и A.A. Коробкина [3], А.Г. Терентьева [4, 5], Э.И. Григолюка и А.Г. Горшкова [6], М.В. Норки-на [7], В.И. Юдовича [8], Г.Г. Шахверди [9], R. Wagner [10], R. Zhao и О. Faltinsen [11], М. Greenhow и W.M. Lin [12], X. Zhu [13], S. Shao [14] и др.
При математическом моделировании подобных задач используется модель сплошной среды, для описания движения которой традиционно применяют два подхода: Эйлера и Лагранжа. Далее рассмотрим численные методы решения задач гидродинамики, основанные на этих подходах, их преимущества и недостатки.
Методы, основанные на подходе Эйлера, используют стационарную, чаще всего регулярную сетку, сквозь которую движутся частицы (малые объемы) сплошной среды, а все физические характеристики определяются в узлах данной сетки, т.с. они не связаны с конкретными материальными частицами, а в каждый момент времени являются характеристиками разных частиц, находящихся в данный момент в данной точке пространства. Методы этого класса позволяют рассчитывать задачи с большими деформациями и чаше всего применяются для задач гидро- и газодинамики. Решения, полученные с их помощью, обладают высокой точностью, кроме того они хорошо изучены
и имеют проработанное теоретическое обоснование. Сложность применения методов данного класса к решению задач со свободными границами является заранее неизвестное положение свободной границы и вытекающие отсюда проблемы с постановкой граничных условий.
Методы, основанные на Лагранжевом подходе, используют подвижную сетку, которая представляет собой дискретное представление материальной среды. Узлы такой сетки жестко связаны ребрами и вместе с ними образуют ее ячейки. Б этом случае сетка двигается и деформируется вместе со сплошной средой, при этом связи узлов сохраняются. Физические характеристики, определяемые в узлах сетки, являются характеристиками соответствующих частиц материальной среды. В отличие от эйлеровых, данные методы позволяют легко отслеживать свободные границы и границы раздела, но также имеют и недостатки, наиболее значительный из которых - невозможность рассчитывать задачи с большими деформациями расчетной области, поскольку они приводят к значительным деформациям расчетной сетки вплоть до пересечений границ (ребер) ячеек, что, в свою очередь, влечет за собой аварийное завершение программы.
Существуют также методы, основанные на совместном использовании обоих подходов, разрабатываемые с целыо устранения вышеупомянутых недостатков.
Самыми распространенными и хорошо изученными изданный момент являются методы конечных разностей (Finite Difference Methods, FDM) [15-18]. Они основаны на эйлеровом подходе, а для получения разностных схем решаемых дифференциальных уравнений используют разложение характеристик, входящих в эти уравнения, в ряды Тейлора. На данный момент существует огромное разнообразие разностных схем, с их помощью решено и до сих пор решается большое количество прикладных задач, хорошо проработана теория, изучены аппроксимационные характеристики схем, их устойчивость и сходимость.
7
К разряду эйлеровых также относится метод конечных (контрольных) объемов (Finite Volume Method, FVM) [19,20], основанный на интегральных законах сохранения. На первом этапе для любого конечного объема формулируется закон сохранения. Затем расчетная область покрывается сеткой, в узлах которой будут рассчитываться физические характеристики моделируемого процесса. Далее выбираются контрольные объемы, чаще всего, с центрами в узлах расчетной сетки и границами, проходящими через центры ребер ячеек сетки. Для каждого полученного контрольного объема записывается дискретный аналог закона сохранения на основе баланса всех потоков через границы рассматриваемого объема. Метод конечных объемов в большинстве случаев позволяет получать консервативные схемы, допускает дискретизацию расчетных областей со сложной геометрией, а также позволяет строить более точные схемы вблизи 1раниц области по сравнению с методами конечных разностей. Эги достоинства метода обусловлены возможностью использовать нерегулярные сетки, равно как и контрольные объемы произвольной формы. Отличительной особенностью данного метода является то, что законы сохранения применяются на этапе построения численных схем, а не на более раннем этапе вывода дифференциальных уравнений, как, например, в методах конечных разностей. Кроме того, физические законы сохранения выполняются не в предельно малых объемах (частицах) среды, а в конкретных конечных подобластях.
Для отслеживания свободной поверхности или контактных границ метод конечных объемов может комбинироваться с методом Volume of Fluid (VOF). Метод VOF был разработан в Национальной лаборатории Лос-Аламоса (США) в конце 70-х - начале 80-х годов [21]. Одна из главных особенностей метода - возможность расчета течений в многосвязных областях с наличием разрывов характеристик и больших деформаций свободной поверхности. В данном методе в качестве маркера, позволяющего определять положение
8
свободной поверхности, служит функция объемной концентрации среды в ячейке [22].
В середине 50-х годов появился метод частиц в ячейках (Particle-in-Cell, PIC), разработанный группой ученых во главе с Харлоу [23] также в лаборатории Лос-Аламоса. Метод сочетает в себе оба подхода к описанию движения сплошной среды - существует как неподвижная эйлерова сетка, так и набор движущихся сквозь нее лагранжевых частиц. На эйлеровом этапе рассчитываются предварительные значения скоростей с учетом лишь вклада давления, затем, на лагранжевом этапе, рассматривается поток частиц через границы ячеек и, таким образом, учитывается вклад конвективных членов [20,23,24]. Давление вычисляется из уравнения состояния. Метод предназначен для моделирования течений сжимаемой среды, однако, в силу постоянства массы частиц, уравнение неразрывности (сохранения массы) во внимание не принимается. Несмотря на то, что изначально метод был разработан для решения уравнений Эйлера, он также позволяет решать уравнения движения при наличии вязкого трения. В зависимости от решаемой задачи в уравнения также можно включить и искусственную вязкость. Кроме того, расчеты можно проводить в любой ортогональной криволинейной системе координат [23]. Метод позволяет рассчитывать многофазные течения без каких-либо ограничений на степень деформации границ раздела и свободных поверхностей.
Несмотря на то, что метод частиц в ячейках позволил значительно расширить класс моделируемых численными методами физических явлений, он не был свободен от недостатков. В частности, результаты расчетов, полученные с его помощью при относительно малом числе расчетных частиц были неточны - наблюдались значительные осцилляции гидродинамических величин. Использование же в расчетах большего числа частиц оказывалось неразумным для ЭВМ того времени. В связи с этим были разработаны экономичные модификации метода, а именно метод жидкости в ячейках (Fluid-in-Cell, FLIC) [25] и метод крупных частиц [26]. Упомянутые методы являются
/
весьма схожими и на эйлеровом этапе не отличаются от метода Харлоу. На лагранжевом же этапе вместо перемещения дискретного набора частиц используется поток массы через границы ячеек. Еще одна модификация - метод Fluid-lmplicit-Particle (FLIP) [27], представляющий собой обобщение метода РТС на случай подвижной адаптивной эйлеровой сетки в целях повышения локальной точности решения. Наиболее полный материал по методам частиц в ячейках, содержащий как математические основы метода, так и конкретные его приложения в современных расчетах, включая коды ирофамм на языке Фортран можно найти в монографии [24]. Хороший обзор метода PIC, его модификаций и приложений содержится также в работе [28].
Для расчета течений несжимаемой жидкости со свободными границами в 1965 году, также иод руководством Харлоу, был создан метод маркеров и ячеек (Marker-and-Cell, MAC) [29]. В отличие от метода Р1С частицы здесь в расчете физических характеристик не участвуют, представляя собой лишь маркеры свободной границы. В методе применяется схема расщепления, а для определения давления необходимо решать уравнение Пуассона. Схема расщепления метода MAC очень схожа с проекционной схемой Чорина [30]. Позднее появились модификации метода, такие как Simplified MAC (SMAC) [31], Stanford-University-Modificd MAC (SUMMAC) [32], Semi-Implicit MAC (SIMAC) [33], MAC-Reynolds-Low (MACRL) [34] и др.
Метод функций уровня (Level Set Method) был предложен Osher и Sethian в 1988 году [35]. В качестве маркера свободной поверхности или границы раздела вместо дискретного набора частиц служит линия уровня некоторой функции. Достоинством метода является простога описания различного рода кривых, прямое вычисление геометрических характеристик свободных поверхностей и границ раздела - кривизны, касательной, нормали и т.д. Кроме того, применение метода к решению трехмерных задач не составляет каких-либо дополнительных трудностей. Подробное описание метода Level Set можно найти в монографиях [36,37].
10
Начиная с 50-х годов прошлого века значительное распространение получил метод конечных элементов (Finite Element Method, FEM) [38-41]. Расчетная область в МКЭ представляет собой сетку, узлы которой сохраняют между собой жесткие связи и двигаются вместе с материальной средой, а ячейки сетки в методе принято называть элементами. Его достоинства заключаются в легком внедрении граничных условий, достаточно высокой точности и возможности проследить всю эволюцию свободной границы. К достоинствам можно также отнести наличие проработанной теоретической базы, большое количество доступной литературы как по теории метода, гак и по его приложениям, а также широкий спектр решаемых им задач. Основной недостаток, как и у других сеточных лагранжевых методов - невозможность проведения расчетов в областях со сложным поведением свободной поверхности в силу перехлеста границ ячеек расчетной сетки за счет сильных деформаций расчетной области. Наиболее полное описание метода и его приложений, как к задачами упругости, так и к задачам гидродинамики, можно найти в 3-х томной монографии Зенкевича [42-44].
Известное распространение также получили метод фаничных элементов (Boundary Element Method, BEM) [45-47] и комплексный метод граничных элементов [48,49]. Удобство использования данных методов обусловлено тем, что консчноэлсментной дискретизации подвергается лишь граница расчетной области, поэтому такие элементы называются граничными. При этом, в любой точке области решение может быть получено по значениям на границе. Данные методы, однако, обладают все тем же недостатком - невозможность расчета задач с сильными деформациями границы расчетной области.
Метод Arbitrary Lagrangian-Eulcrian (ALE) [50], как и ранее рассмотренные методы типа PIC и MAC, относится к комбинированным лафанжево-эйлеровым методам, что напрямую следует из его названия. В данном методе используется сетка, которая может двигаться произвольно, т.е. она не остается фиксированной, как эйлерова сетка, но и не подчиняется законам
11
движения лагранжсвой сетки, как, например, в методе конечных элементов, откуда и слово "произвольный" (arbitrary) в названии метода. Алгоритм метода разбивается на три основных этапа [51]: 1 - перемещение сетки; 2 -перестройка сетки; 3 - интерполяция значений со старой сетки на новую. Удобство данного метода заключается в том, что можно перестраивать сетку лишь в тех местах, где это необходимо, например, там, где лагранжева сетка сильно деформирована, что в свою очередь могло бы повлиять на точность получаемых результатов, либо вообще способствовать аварийному завершению работы алгоритма. Недостатком метода является наличие процедуры интерполяции, которая способствует сглаживанию результатов. Кроме того, для отслеживания свободных поверхностей и границ раздела типа "жидкость-жидкость "жидкость-тело сетка в методе ALE вблизи них должна вести себя подобно обычной лагранжевой сетке, что в свою очередь приводит к общему недостатку сеточных лагранжевых методов - перехлесту границ элементов. Данный метод часто используется для решения задач о взаимодействии жидкости с погруженными телами.
Обозначенные недостатки сеточных методов, в особенности методов лагранжевой природы, привели к появлению так называемых бессеточных методов, которые в последнее время получают все большее распространение. Несмотря на употребляемый термин, слсдуег, тем не менее, понимать, что не все методы, относящиеся к классу бессеточных, вообще не используют сетку при расчетах. Поскольку строгого, устоявшегося определения бессеточных методов нет, будем следовать терминологии Лью [52] и определять бессеточные методы как методы, не требующие использования связной сетки, по крайней мере, для построения функций формы (например, в методе конечных элементов для данной процедуры обязательно используется сетка). К таковым относятся бессеточные методы, основанные на слабой форме уравнений, поскольку для ее интегрирования все же требуется наличие сетки. Идеальным же требованием к бессеточному методу является отказ от сетки на любом
этапе численного решения задачи. К данному классу относятся методы, использующие дифференциальную форму уравнений механики жидкости.
Основное отличие бессеточных методов от классических лагранжевых состоит в том, что сетка строится на каждом временном шаге по новому набору узлов. Это означает, что, по ходу проведения вычислений, узлы области расчета могут перемещаться свободно, в силу отсутствия между ними жестких топологических связей. Такой подход приводит к ряду преимуществ бессеточных методов перед традиционными сеточными при решении задач с большими деформациями расчетных областей:
- отсутствие необходимости в применении сложных и ресурсоемких алгоритмов адаптации сетки, с целью избежания самопересечения ее ребер, что в обычных случаях приводит к аварийному завершению расчетов, а также необходимости использования стандартной в таких случаях процедуры интерполяции, результатом которой является неизбежное понижение точности результатов.
- возможность рассчитывать задачи с разрывами характеристик с заведомо более высокой точностью, в силу того, что поверхность разрыва не должна проходить строго по границам элементов, как того требуют методы конечных элементов.
- возможность использования простых адаптивных процедур добавления и удаления узлов в локальных областях, поскольку в бессеточных методах такая процедура не влечет за собой нарушения связности сетки.
Далее более подробно рассмотрим основные особенности наиболее распространенных численных методов, относящихся к классу бессеточных.
Метод Moving Least Squares (MLS) был предложен в работе Lancaster и Salkauskas [53]. Аппроксимация функции в точке находится путем минимизации функционала, представляющего собой сумму взвешенных квадратов отклонений значений аппроксимированной функции от точных значений этой
13
•г
функции в узлах сетки. В классическом методе наименьших квадратов в качестве весовой функции используется функция-константа равная 1, в методе MLS - функция, имеющая форму гауссовой кривой и обладающая компактным носителем. Глобальная аппроксимация функции получается путем перемещения точки максимума весовой функции по точкам аппроксимации в пределах расчетной области. Эта особенность метода добавила к названию метода наименьших квадратов слово "перемещение" (moving). Компактность носителя весовой функции позволяет использовать для аппроксимации функции в точке лишь ближайшие узлы, а ее форма дает возможность по-разному оценивать вклад в аппроксимацию различных узлов. MLS не является методом численного моделирования физических процессов, а лишь способом построения функций формы, весьма распространенным во многих бессеточных методах, в частности, MLSPIl, DEM, EFG, которые будут рассмотрены ниже. Основной недостаток заключается в том, что функции формы, построенные на основе метода MLS, не удовлетворяют условию Кронекера (см.
2.3.), что является причиной общей для всех методов, построенных на MLS, проблемы - сложности внедрения главных граничных условий, поскольку значения аппроксимированной функции не совпадают с точными значениями даже в узловых точках. Существует также метод MLS для комплексных переменных [54].
Разумеется, метод MLS, несмотря на довольно широкое его использование при построении функций формы в бессеточных методах, не является единственным инструментом для этих целей. Liu и Gu [55] предложили новый численный метод - метод точечной интерполяции (Point Interpolation Method, PIM), использующий слабую форму уравнений, процедура построения функций формы в котором отличается от MLS. Новая процедура устраняет принципиальный недостаток метода MLS - нарушение условия Кронекера. Другим отличием является то, что размер полиномиального базиса (количество одночленов в базисе) для метода PIM должен соответствовать количеству точек в
14
области-носителе, на которой строится функция формы. К общим недостатком процедур MLS и PIM относится необходимость обращения, во многих случаях, вырожденных матриц.
Как уже было сказано ранее, бессеточные методы, основанные на слабой форме уравнений, требуют для целей интегрирования построения глобальной сетки на всей расчетной области, что является существенным недостатком этого класса методов. С целью ею устранения метод точечной интерполяции был модифицирован для локальной слабой формы уравнений [56], что подразумевает построение сетки интегрирования лишь в малой подобласти расчетной области. Разумеется, этот процесс является значительно менее трудоемким. Далее будут рассмотрены другие методы, также основанные на локальной слабой форме уравнений.
Несмотря на удобство использования полиномиальных базисов в методе PIM, матрица, полученная в ходе построения интерполяции и требующая последующего обращения, может оказаться вырожденной, что является принципиальной проблемой, не позволяющей продолжить процесс вычислений. В связи с этим, было предложено [57] вместо полиномиальных базисных функций использовать радиальные, т.е. зависящие не от координат узлов, а лишь от расстояния между ними. Использование радиальных базисных функций приводит к симметричной матрице, всегда обратимой. Модификация метода PIM, использующая радиальные базисные функции получила название RPIM (Radial PIM).
Ранее было сказано, что методы, построенные на базе граничных интегральных уравнений (МГЭ, КМГЭ), имеют ряд преимуществ перед методами конечных элементов. Идея использования граничных интегральных уравнений для построения новых методов, по аналогии с методами граничных элементов, привела к появлению метода граничной точечной интерполяции (Boundary PIM) [58]. Как следует из названия, для построения функций
15
- Киев+380960830922