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

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

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

Лабораторное моделирование и исследование поверхностных волн аналитическими методами широко используются для проверки достоверности численных результатов при проведении вычислительных экспериментов, для определения областей применимости тех или иных математических моделей, для тестирования численных методик и для получения начальной информации об изучаемом явлении. В работах [9, 16, 17, 45] приведены экспериментальные данные о влиянии параметров волн на величину максимального заплеска на вертикальную преграду, рассмотрен накат волн на некоторые гидротехнические сооружения, дана оценка воздействия волн на препятствия. Авторами работ [127, 139] проведены эксперименты о прохождении уединенных волн над различными донными препятствиями. Обзор экспериментальных работ, посвященных исследованию поверхностных волн, приводится в работе [150].

Аналитические решения задач о построении форм уединенных волн, их взаимодействии с препятствиями и трансформации при распространении приведены в работах [40, 42, 49, 63, 65, 102, 122]. В монографиях [38, 42, 63, 102] излагается теория волновых движений жидкости и ее приложение, содержится разбор специальных вопросов этой теории, относящихся к ряду задач геофизики, теории движения кораблей и др. Рассматриваются различные нелинейные эффекты, возникающие при распространении и взаимодействии волн с погруженными телами. В [102] содержится обзор развития теории нелинейных волн за три последних десятилетия.

В работе [122] исследуется приближение соленоидальных и кноидальных волн на основе теории мелкой воды. В приведенной модели отсутствует предположение о бесконечно малом перемещении жидкости по вертикали, которое делалось в более ранних работах других авторов. Построение уединенных волн также является одной из центральных тем работы [49]. В работах [40, 65] рассматривается движение жидкости со свободной поверхностью. Исследовано косое набегание волн малой амплитуды на упругую пластину, плавающую на поверхности жидкости. Определены коэффициенты отражения и прохождения волн, а также вертикальные смещения пластины.

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

В задачах, имеющих практический интерес, таких как: движение волн по поверхности, возникновение цунами, движение тел в жидкости, размеры длин волн обычно бывают большими, поэтому при исследовании распространения волн и их взаимодействия с твердыми преградами, можно пренебрегать влиянием вязкости и рассматривать задачи в постановке, основанной на модели идеальной несжимаемой жидкости [44, 60, 61, 74].

Выбор математической модели для исследования поверхностных волн зависит от особенностей течения. Достаточно большая часть используемых в настоящее время численных методов расчета поверхностных волн основана на применении тех или иных приближенных математических моделей, в качестве которых обычно берутся различные приближения теории мелкой воды. В этих моделях, как правило, пренебрегают изменением параметров жидкости по глубине, что может дать в некоторых случаях, например, вблизи твердой стенки, некоторые погрешности. В рамках данной модели решен класс задач, описанный в [23, 47, 70, 74, 77, 113, 140].

Численное моделирование уединенных волн на основе дискретной модели несжимаемой жидкости описывается в работах [70, 72]. Приводятся результаты расчетов уединенных волн различной амплитуды, кинематические и динамические характеристики наката волн на вертикальную и наклонные стенки. В работе [113] исследуется распространение волн по поверхности жидкости в бассейне переменной глубины. Авторами работ [47, 77] разработан инструментарий вычислительного эксперимента для решения задач проблемы цунами, включающий систему математических моделей и алгоритмов для численного моделирования распространения и трансформации волн цунами и определения времени добегания волн цунами от очаговой зоны до побережья. Решены модельные и прикладные задачи проблемы цунами и выявлены фундаментальные характеристики, определяющие процессы распространения и трансформации таких волн [77, 140].

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

Большое распространение получили модели потенциальной жидкости в полной нелинейной постановке для решения задач с сильно нелинейными режимами движения волн. В работах [2, 58, 64, 67, 84, 85, 109, 135, 137] моделировались процессы генерации волн, трансформации волн при их движении, взаимодействии волн с различными преградами. Работы авторов [2, 64, 67] посвящены решению стационарных и нестационарных задач со свободными границами. Рассматриваются задачи о движении твердых тел в жидкости со свободными границами, моделировании и распространении нелинейных уединенных волн по бассейну с ровным дном и над дном с препятствием. Проводится анализ волновых движений, возникающих на поверхности невязкой несжимаемой весомой жидкости при взаимодействии волн с вертикальными и наклонными неподвижными границами. В работах.

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

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

Автором работы [137] предлагается эффективный численный метод на основе смешанной Лагранжево-Эйлеровой постановки для моделирования движения нелинейных волн до момента опрокидывания гребня волны. В работе [109] предложен численный метод на основе рядов Фурье для решения задач в полной нелинейной постановке о движении уединенных необрушающихся волн и взаимодействии волн между собой. В работе [58] исследуется трансформация уединенной волны при ее движении над подводным уступом.

Переход от модели потенциального движения к уравнениям Эйлера [44, 60, 61] расширяет круг задач, доступных для численного исследования. В частности, можно рассчитывать и вихревые движения жидкости со свободной поверхностью, обеспечивая завихренность либо надлежащими начальными условиями, либо непотенциальной массовой силой. К тому же систему уравнений Эйлера можно рассматривать как промежуточный этап на пути от потенциальной модели к уравнениям Навье-Стокса. Обзор методов решения различных задач идеальной жидкости со свободными границами можно найти в монографиях [38, 41, 63].

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

— при помощи аппроксимации дифференциальных операторов в уравнениях более простыми локализованными алгебраическими операторами, действующими в последовательностях узлов, находящихся в области — метод конечных разностей (МКР) [26];

— при помощи представления самой области элементами среды, не являющимися бесконечно малыми, которые в совокупности аппроксимируют реальную систему — метод конечных элементов (МКЭ) [24, 39] и метод контрольных объемов (МКО) [26, 51].

Также большое распространение получили методы граничных элементов (МГЭ) [13, 15] и комплексный метод граничных элементов (КМГЭ) [20]. Привлекательность этих методов обусловлена, прежде всего, тем, что в МКР, как и МКЭ, требуется разбиение всей области течения, в то время как в методе граничных элементов дискретизации подвергается лишь граница расчетной области. Для реализации такой возможности в МГЭ и КМГЭ требуется переход от исходной краевой задачи для дифференциальных уравнений, описывающих некоторый процесс, к соотношениям, связывающим неизвестные функции на границе области. МГЭ и КМГЭ широко используются для моделирования движения поверхностных волн, но они пригодны для расчета только потенциальных течений.

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

Также большое распространение получила группа методов, берущая начало от метода Харлоу — метода частиц в ячейках [75]. Позже появились различные модификации этих методов — методы крупных частиц [12] и метод частиц-в-ячейках [19, 81]. Это направление методов сочетает в себе в определенных чертах преимущества лагранжевого и эйлерового подходов. Область решения разбивается неподвижной (эйлеровой) сеткой, а сама сплошная среда представляется совокупностью частиц (лагранжева сетка), которые движутся через эйлерову сетку ячеек. Частицы служат для определения параметров самой жидкости (энергии, скорости), в то время как эйлерова сетка используется для определения параметров полей давления, плотности, температуры. Главный недостаток этой группы методов заключается в относительно невысокой точности, при которой обычный уровень погрешности составляет несколько процентов.

С ростом производительности компьютеров возродился интерес к лагранжевому описанию среды на основе свободно-лагранжевых методов. В таких методах отношение соседства частиц не фиксировано и может со временем изменяться. Данные методы известны как бессеточные методы [88, 112, 124, 125], которые аппроксимируют уравнения в частных производных, основываясь только на наборе узлов без использования расчетной сетки. Характерными представителями этой группы методов являются метод сглаженных частиц (Smoothed Particle Hydrodynamics, SPH) [131, 132], полунеявный метод движущихся частиц (Moving Particle Semi-implicit, MPS) [120, 121], метод лагранжево-эйлеровых частиц [71].

Метод SPH был предложен в 1977 году Льюси и Монаганом для решения задач астрофизики, позднее метод был применен для решения задач гидродинамики [131, 132]. Основная идея метода состоит в дискретизации среды конечным набором лагранжевых частиц, которые движутся со скоростью потока и допускают произвольную связность между собой. Все функции, входящие в систему дифференциальных уравнений, представляются в виде интегралов по области течения с весовой функцией Дирака. Далее функция Дирака заменяется финитной функцией ядра W, которая может быть продифференцирована аналитически, а интегралы заменяются конечной суммой по соседним частицам. Так как функция ядра легко дифференцируема, то вычисление производной от искомой функции не составляет труда.

Метод MPS [120, 121] является одной из модификаций метода SPH. Основные идеи метода схожи с идеями метода сглаженных частиц, однако он предназначен для расчета течений несжимаемой жидкости. Функция ядра W должна удовлетворять всем свойствам, что и в методе SPH, кроме требования j* WdQ = 1, и зависит от количественной плотности частицы, с помощью п которой достигается условие несжимаемости жидкости. Из-за специфического вида функции ядра производная от искомой функции в методе MPS вычисляется по формуле производной по направлению [120].

Метод лагранжево-эйлеровых частиц впервые описан Франком A.M. в 1992 году [71]. Рассматриваемый метод частиц базируется на свободно-лагранжевой модели, которая строится на основе принципа наименьшего принуждения Гаусса. С другой стороны, массовые, внутренние и поверхностные силы вычисляются в эйлеровых координатах. Полученная в результате численная схема является полностью консервативной и устойчивой. Метод представляет собой некоторый специальный вариант метода Галеркина, и основная идея заключается в переходе от дифференциальной постановки задачи к вариационной. Сложность данного метода заключается в выборе базисных функций, с помощью которых определяются неизвестные функции.

К общим недостаткам перечисленных бессеточных методов можно отнести сравнительно невысокую точность и трудность введения граничных условий [112, 130]. Для внедрения граничных условий используются различные методы: метод множителей Лагранжа, вариационные методы, метод штрафов. В некоторых бессеточных методах для внедрения граничных условий строится специальный вид функций формы. Также большое распространение получил способ внедрения граничных условий на основе разбиения области конечными элементами и применения стандартной процедуры МКЭ, изменяющей коэффициенты матрицы. Дополнительно в бессеточных методах для аппроксимации неизвестных функций требуется обеспечить узловую связность и вычислить значение производной от функции формы, что в некоторых методах представляет серьезные трудности.

В связи с вышесказанным, особенно актуальным является разработка вычислительных алгоритмов, объединяющих сеточный и бессеточный подходы и использующие преимущества каждой из методологий. Характерными представителями этой группы методов являются бессеточный метод конечных элементов (Meshless Finite Element Method, MFEM) [103, 133] и метод естественных соседей (Natural Element Method, NEM) [134].

Методы NEM и MFEM основываются на интерполяции естественных соседей (natural neighbour interpolation) [11, 87, 141, 142]. В вычислительной геометрии понятие естественных соседей связано с понятием ячейки Вороного: для заданного узла pj ячейка Вороного — объединение точек пространства, расстояние от которых до pt меньше, чем до остальных узлов [52, 90]. Ячейки Вороного разбивают все пространство на выпуклые многоугольники, в каждом из которых находится только один из узлов системы точек. Среди многочисленных свойств этих ячеек особо следует отметить единственность разбиения, двойственность триангуляции Делоне и непрерывную зависимость всех геометрических параметров ячеек от координат узловых точек. Естественными соседями точки р{ являются такие точки, ячейки Вороного которых имеют общее ребро. На основе диаграммы Вороного можно разработать эффективные алгоритмы поиска естественных соседей, которые применяются при вычислении функций формы для методов NEM и MFEM.

В основе MFEM лежит классический метод конечных элементов. Для аппроксимации неизвестных функций используются функции формы Лапласа [11, 87]. Расчетная область дискретизируется элементами расширенной триангуляции Делоне, обеспечивающей единственность разбиения и, тем самым, единственность результата интерполяции.

Метод естественных соседей был предложен в 1994 году в работе [147] для решения задач теории упругости. В основе метода лежит интерполяция Сибсона, которая строится по произвольно заданному множеству расчетных узлов на ячейках Вороного первого и второго порядков. Коэффициенты интерполяции Сибсона в плоском случае определяются как отношение площади пересечения ячеек Вороного второго порядка для точки х, введенной в заданное разбиение, с площадью ячейки Вороного первого порядка, содержащую точку х. NEM представляет собой разновидность метода Галеркина, в котором базисные и весовые функции совпадают. Для формирования дискретной системы уравнений используется метод взвешенных невязок в слабой форме. Главным преимуществом интерполяции Сибсона является ее четкая определенность, независимость построения и устойчивость на неравномерном распределении расчетных узлов. Носителем функции формы Сибсона является объединение всех кругов Делоне для заданного узла pi.

Коэффициенты интерполяции Сибсона соответствуют размерности рассматриваемого пространства (расстояние в R1, площадь в R2, объем в.

R3), в то время как функции формы Лапласа рассматривают меры, на единицу меньшей размерности. В работах [11, 87] показано, что интерполяция Сибсона позволяет получить более точные результаты, чем интерполяция Лапласа. Однако построение сибсоновских функций формы является трудоемкой задачей вычислительной геометрии, что увеличивает временные затраты при вычислении коэффициентов интерполяции. В работах [92, 149] предложен алгоритм Бовье-Вотсона построения интерполяции естественных соседей, который основывается на разбиении области ориентированными треугольниками и минимизирует время вычисления коэффициентов. Для внедрения граничных условий используется подход, который применяется в МКЭ. Тем самым в методе естественных соседей преодолеваются основные трудности, свойственные большинству бессеточных методов, — вычисление функций формы, вычисление производных от искомой функции и внедрение граничных условий.

Многие известные бессеточные численные методы (SPH, MPS, PIM, VOF) [124, 125] успешно применяются для моделирования задач гидродинамики с большими деформациями расчетной области. Данные методы позволяют достаточно точно воспроизводить кинематику течений, однако получение динамических характеристик, необходимых для расчета гидродинамических нагрузок, является весьма затруднительной задачей. Актуальность определения динамических нагрузок при взаимодействии поверхностных волн с береговыми и донными сооружениями обусловлена вероятностью появления катастрофических последствий, возникающих в том случае, когда гидродинамические нагрузки превышают допустимые пределы. При проектировании морских сооружений необходимо оценивать степень воздействия поверхностных волн на различные объекты. Проблема прочности конструкций и сооружений, взаимодействующих с жидкостью, занимает центральное место при оценке их эффективности и времени жизни.

Благодаря своей точности и устойчивости NEM нашел приложение в решении дифференциальных уравнений [93], в механике твердого тела [143, 144], в теории упругости [96, 101, 145]. При решении данного класса задач деформации расчетной области незначительны, поэтому, несмотря на все свои преимущества, применение метода NEM к решению задач механики жидкости со свободными границами затруднительно. В данной работе проводится модификация метода естественных соседей для решения задач гидродинамики с большими деформациями расчетной области. Разработка численных методов для моделирования течений идеальной несжимаемой жидкости, позволяющих с высокой точностью определять не только кинематические, но и динамические характеристики течений с большими деформациями расчетной области, является важной и актуальной задачей современной гидродинамики.

О предмете и содержании диссертации.

Диссертация состоит из введения, четырех глав, заключения и списка литературы. Главы диссертации разбиты на параграфы, нумерация которых внутри каждой главы своя. Общий объем диссертации составляет 163 листа.

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

состоит из 150 наименований.

Основные результаты диссертации состоят в следующем:

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

2. Разработан и реализован в виде программного комплекса численный алгоритм модифицированного метода естественных соседей для решения в полной нелинейной постановке плоских нестационарных задач со свободными границами. В качестве схемы движения по времени реализован метод дробных шагов. Для решения построенной системы линейных алгебраических уравнений реализован метод сопряженных градиентов с 01Ьи (0)-предобусловливанием, который значительно сокращает временные затраты на проведение численных экспериментов.

3. Реализован алгоритм «заметающей плоскости» представления расчетной области ячейками Вороного. Разработаны на основе дискретизации области ячейками Вороного алгоритмы определения естественных соседей для узловой точки и формирования структуры межузловой связности, а так же алгоритм определения границ многосвязной расчетной области.

4. Разработан и реализован в виде программного комплекса параллельный алгоритм метода естественных соседей для многопроцессорных вычислительных систем с распределенной памятью.

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

Заключение

.