Содержание к диссертации
Введение
Глава 1. Общая постановка задачи и метод решения 16
1 Постановка задачи 16
2 Метод движения по времени 19
3 Построение сетки граничных элементов 22
3.1 Первый подход 23
3.2 Второй подход 26
4 Метод граничных элементов 28
4.1 Граничные элементы и интерполирующие функции ... 30
4.2 Вычисление интегралов 33
4.3 Вычисление диагонального коэффициента матрицы Н. 39
4.4 Функция Грина для полубесконечных областей 42
4.5 Решение системы линейных алгебраических уравнений 44
5 Вычисление поля скоростей 46
6 Задача о движении абсолютно твердой сферы в
безграничной идеальной несжимаемой жидкости 48
Глава 2. Динамика сферического пузыря 52
1 Задача Релея 52
2 Замыкание газонаполненного пузыря 57
3 Закон сохранения энергии 61
Выводы 66
Глава 3. Влияние силы тяжести, положения твердой стенки и поверхностного натяжения на динамику пузыря 67
1 Вычисление характеристик струйки. Метод оценки урона, наносимого струйкой твердой стенке 67
2 Влияние силы тяжести 72
3 Влияние твердой стенки 80
3.1 Плоские стенки 80
3.2 Угловые стенки 87
4 Влияние поверхностного натяжения 90
5 Сравнительный анализ расчета с использованием пространственной и осесимметричной моделей 93
6 Сравнительный анализ расчетов полученных с помощью разработанного метода граничных элементов с расчетами описанными в работе Вэнга (Q.X. Wang) 96
Выводы 101
Глава 4. Динамика парогазовых пузырей около наклонных стенок 103
1 Импульс количества движения 103
2 Исследование форм пузыря и характеристик струйки 108
2.1 Угол наклона стенки є = л 708
2.2 Угол наклона стенки с=Зж/4 116
2.3 Угол наклона стенки є=я!і 123
2.4 Угол наклона стенки є=яі4 130
2.5 Угол наклона стенки - = 0.0 737
3 Характеристики струйки. Переход к размерным величинам 144
Выводы 148
Литература 150
- Граничные элементы и интерполирующие функции
- Замыкание газонаполненного пузыря
- Сравнительный анализ расчета с использованием пространственной и осесимметричной моделей
- Характеристики струйки. Переход к размерным величинам
Введение к работе
Работа посвящена исследованию динамики пространственного парогазового пузыря в идеальной несжимаемой жидкости. В ней рассматривается эволюция паровой и парогазовой сферической каверны в безграничной идеальной несжимаемой жидкости, изучается влияние силы тяжести, газосодержания и поверхностного натяжения на процесс эволюции пузыря. Основное внимание уделяется изучению взаимодействия парогазовой полости с различными твердыми стенками. В качестве инструмента численного исследования используется метод граничных элементов.
Описанная в данной работе модель применяется для моделирования различных, непохожих на первый взгляд, явлений. С одной стороны рассматривается динамика подводных взрывов, с другой стороны прогнозируется урон от кавитациошюй эрозии, когда изучение всего явления пузырьковой кавитации невозможно без изучения динамики одной, отдельно взятой, кавитационной полости. В последнее время особую актуальность приобрели задачи биофизики, а именно задача об исследовании течений в кровеносных сосудах [20, 21], а также задачи, возникающие при изучении надежности работы протеза митрального клапана сердца [86], в этом случае оценивается возможность урона от пузырьков жидкости, находящихся в потоке крови, стенкам митрального клапана. Что же объединяет эти столь разные на первый взгляд явления? Одна из таких черт - это возникающие кумулятивные эффекты. Пузырь, развиваясь из кавитиционого зародыша (если это кавитационный пузырек), или образующийся при взрыве заряда (если рассматривается подводный взрыв), в процессе своего роста, как правило, сохраняет форму близкую к сферической [52], Достигнув наибольшего объема, пузырь переходит в фазу замыкания. Близость твердой границы и (или) действие силы тяжести нарушают одномерность течения, даже если в момент максимального расширения полость была сферической. В ряде случаев в процессе замыкания пузыря формируется струйка жидкости, внедряющаяся в пузырь до момента касания его противоположной стенки. Такая струя может быть направлена в сторону стенки и иметь скорость порядка сотен, а при особых условиях, даже тысяч метров в секунду [102]. Частицы на дальней от стенки поверхности пузыря, в случае если струя направлена к стенке, получают большее ускорение, то есть возникает
классический кумулятивный эффект [47]. Предполагается, что механизм разрушения мишеней определяется именно воздействием высокоскоростной кумулятивной струи, формирующейся на стадии замыкания пузыря [47]. Поэтому актуальной задачей является выяснение роли различных факторов, таких как наличие газа в пузыре, действие силы тяжести, положение твердой стенки, в процессе образования и развития кумулятивной струи.
Исследование задачи об эволюции пузыря имеет достаточно долгую историю. Особый вклад в понимание механизма развития сферических каверн, представляющих собой идеализированную модель нестационарных перемещающихся каверн, внес лорд Релей, опубликовавший в 1917 г. статью "О давлении развивающемся в жидкости при схлопывании сферической каверны". Релей использовал предложенную в 1859 г. Безантом постановку задачи о сферически симметричной пустой полости, мгновенно образующейся в бесконечном объеме невесомой идеальной несжимаемой жидкости, при постоянном давлении на бесконечности. Теоретическое исследование стационарного безотрывного обтекания сферического пузыря в вязкой безграничной жидкости нашли свое отражение в работах [26, 60] и др.. Результаты этих исследований показали, что существует широкий диапазон чисел Рейнольдса, при которых распределение нормальных напряжений на границе пузыря мало отличается от давления в жидкости вблизи его границы при потенциальном обтекании потоком идеальной жидкости. Поэтому, исследование стационарной формы пузыря и ее устойчивости при движении можно проводить используя модель идеальной жидкости. В такой постановке задачи о движении пузыря решаются начиная с работы Н.Е.Жуковского [42]. Подробный обзор теоретических и экспериментальных результатов, полученных до 1975 года, приведены в работах Е.И.Забабахина [43], а также О.В.Воинова и А.Г.Петрова [25]. Достаточно полно рассматривается динамика сферических кавитационных пузырей в книге Ю.Л.Левковского [54], влияние твердой стенки на замыкание сферической каверны рассмотрено в работе ЮЛ.Левковского и Г.Г. Судакова [55]. Последние достижения в этой области отражены в работе В.К.Кедринского [47]. Поведение пузыря в вязкой жидкости описывается в работе П.К.Волкова [28].
Решение нестационарных задач о движении и схлопывании пузыря в жидкости приводит к необходимости рассматривать нелинейные эффекты, которые характеризуются накоплением кинетической энергии в довольно малом объеме жидкости и резким ростом скоростей на свободной границе. В
процессе замыкания пузырь уменьшается до столь малого объема, что
влияние сжимаемости на процесс замыкания пузыря становится весьма
существенным. Тем не менее, экспериментальные исследования
показывают, что около 90% времени своего существования пузырь находится в состоянии, когда его радиус не достигает той критической величины, хотя и успевает потерять сферическою форму [48, 52].
Задача о динамике пузыря относится к классу задач со свободными границами, которые бурно развиваются последние годы, что связано с огромным прикладным значением достижений в этой области. Однако аналитические исследования подобных задач весьма затруднительны, так как даже в рамках модели идеальной жидкости наличие свободной границы приводит к существенно нелинейным результатам. Кроме того, на свободной границе, как правило, задаются нелинейные краевые условия. Постановки таких задач и методы их решения описаны в книге М.А.Лаврентьева и Б.В.Шабата [53], а также в работах [17, 26, 40, 41, 68] и многих других.
С появлением быстродействующих и высокопроизводительных компьютеров стало развиваться направление численного моделирования. Первые работы, посвященные исследованию эволюции пузыря около твердой стенки, были выполнены в осесимметричной постановке. Оригинальные численные результаты о процессе схлопывания пузырька в безграничной жидкости около твердой стенки получили M.S. Plesset и R.B. Chapman [99], О.В. Воинов и В.В. Воинов [24], M.S. Plesset и A. Prosperetti [100], Т.М. Mitchell и F.G.Hammit [98], Поведение полусферического пузырька, присоединенного своим основанием к твердой стенке изучали А. Shiraa и К. Nakajima [107]. Работы J.R. Blake и др. [82, 83] посвящены исследованию влияния твердой стенки и свободной поверхности на эволюцию первоначально расширяющегося пузыря. Направление миграции пузыря и направление развития кумулятивной струйки в упомянутых работах определяется через импульс количества движения (импульс Кельвина). Теоретические исследования импульса Кельвина можно найти в работах J.P.Best и J.R.Blake [76, 78]. Набор методов, используемых для изучения эволюции как кавитационных, так получаемых в результате взрывов, пузырей, достаточно широк и по сей день постоянно пополняется. Один из методов широко применяемых для данного круга задач - метод граничных элементов, использованный и в данной работе. Наиболее известные работы посвященные моделированию осесимметричного пузыря методом граничных
элементов - [77, 82, 83] и др. Известно всего несколько подобных работ в пространственной постановке [85, 111]. В работе [112] описан подход, позволяющий проводить расчет перехода от пузыря к тору, то есть от односвязной к двусвязной области методом граничных элементов. Использование метода конечных элементов для задач эволюции пузыря описано в работах [103, 109], метода объемов жидкости (volume-of-fluid method) - [97], метода Лагранжа-Томсона - [ПО], метода обобщенных вихрей (generalized vortex method)- [101].
Большое внимание при изучении динамики пузыря отводится также эксперименту, хотя техника эксперимента в этом случае чрезвычайно сложна [79, 80]. Экспериментальные результаты можно найти в работах Benjamin Т.В. и Ellis А.Т. [75], Lauterborn W. и Bolle Н. [96], Shima А. и Nakajiama К. [107], Гривни на Ю.А. Зубрилов СП. [30] а также в ряде более поздних работ Kedrinskii К.А. и Stepanov V.A. [94], Tomita Y. и др.[105], Shima А. [106], Blake J.R. и др. [84], S.J. Shaw [108], N.J. Lawson и др. [97], Field J.E. [87], Hooton M.C. и Blake J.R. [92].
В настоящее время расчет пространственных задач уже не лежит за пределами возможностей современных компьютеров и, более того, не является чем-то экстраординарным. Тем не менее, как сам расчет пространственной задачи, так и анализ результатов все еще остается достаточно сложным и трудоемким процессом, известно всего несколько работ посвященных исследованию эволюции пространственных пузырей [85], взаимодействию пузыря с твердой стенкой [111]. Работа [84] посвящена как экспериментальному так и численному моделированию эволюции пузыря и тора образующегося в результате замыкания пузыря. В работе [101] описывается исследование пульсации пространственного пузыря в безграничной несжимаемой жидкости с учетом поверхностного натяжения в поле переменного давления.
Основная цель данной работы - изучить поведение пространственного парогазового пузыря при его эволюции как в безграничной жидкости, так и вблизи твердой стенки, а также выяснить влияние различных факторов, таких как наличие газа в пузыре, действие силы тяжести, поверхностного натяжения, положения стенки, на процесс образования и развития кумулятивной струи. Кроме того, необходимо оценить урон, наносимый струйкой твердой стенке (мишени), получить размерные характеристики струйки для различных типов пузырей
(кавитационных пузырей и пузырей, возникающих при взрыве различных зарядов).
В первой главе диссертации изложена общая постановка нелинейной пространственной задачи о динамике пространственного парогазового пузыря в идеальной несжимаемой тяжелой жидкости с учетом поверхностного натяжения, описана методика движения по времени и реализация метода граничных элементов (МГЭ) на каждом временном шаге.
Описание общей нелинейной постановки задачи приведено в работе Л.В. Овсянникова [61], обезразмеривание переменных при переходе к численной реализации - в работах [54, 82].
Используемый в работе метод граничных элементов является разновидностью метода граничных интегральных уравнений, который применяется для анализа краевых задач еще с начала века. Неоспоримым преимуществом метода является то, что для нахождения неизвестных функций внутри области используются уравнения, записанные по границе, и являющиеся формулировкой исходной поставленной задачи, ведущей к точному решению. Известно, что численное интегрирование всегда представляет собой более устойчивый и точный процесс, чем численное дифференцирование. Кроме того, переход к уравнению, записанному по границе, позволяет снизить размерность задачи на единицу, что в пространственном случае чрезвычайно важно. Матрица порожденная МГЭ является полнозаполненной, в отличии от других популярных методов — метода конечных элементов [3, 44, 49, 63] и разностных схем [16, 59], при реализации которых матрица получается сильно разреженной или ленточной. Основные теоретические результаты и современные достижения в области решения интегральных уравнений отражены в работах В.Г.Мазьи [58], А.Ф.Верлаия и В.С.Сизикова [23].
Численная реализация метода граничных интегральных уравнений начинается с метода Крылова-Боголюбова, описанного в 1929 году. Изложение этого подхода, который по принятой на сегодняшний день классификации можно отнести к методам с постоянными граничными элементами или методу "панелей", можно найти в работе Л.В.Канторовича и В.И.Крылова [45].
В данном случае в качестве основного соотношения метода
граничных элементов используется третья формула Грина, записанная для
д<р потенциала поля скоростей (р и его нормальной производной q = -—
С(х}р(х) + $д*(х,ї)ф,ї)<П:(ї)= fr*(x,?)q{x,?)dr(4), (уЛ)
Г г
где ^"-фундаментальное решение уравнения Лапласа, q* =—?— - его
нормальная производная, С = со\2к, со - телесный угол, под которым видна
поверхность из точки х.
Граница области Г аппроксимируется набором плоских или криволинейных треугольных или четырехугольных элементов, на которых искомые величины могут быть аппроксимированы постоянными, линейными или квадратичными функциями. В зависимости от этого элементы называются "постоянными", "линейными" или "квадратичными", что в свою очередь определяет число узлов на элементе. Записав уравнение (v.l) для всех узлов границы (называемых точками наблюдения) дискретный аналог основного соотношения МГЭ (v.l) и вычислив необходимые интегральные коэффициенты, получим систему линейных алгебраических уравнений для неизвестных значений искомых функций в узлах границы. Сведение интегральной формулировки исходной задачи к системе линейных алгебраических уравнений подробно описано в работе [65].
Полный обзор метода граничных элементов можно найти в работах [19] и [18].
Для дискретизации интегрального соотношения (v.l) в диссертации используются линейные треугольные элементы. Соотношение (v.l) на каждом шаге по времени приводится к интегральному уравнению Фредгольма первого рода относительно неизвестных значений нормальной производной, которое является некорректным по Адамару [23]. Однако известно, что если интегральные коэффициенты имеют сильную особенность, которая вычисляется в смысле главного значения по Коши, как в нашем случае, то результирующая система линейных алгебраических уравнений имеет диагональное преобладание, при решении которой наблюдается эффект "саморегуляризации" [58].
Для верификации работы МГЭ на каждом временном шаге проводятся тестовые расчеты по движению твердой сферы в безграничном потоке идеальной несжимаемой жидкости [67],
Вторая глава посвящена решению задачи динамики пузыря в безграничном объеме идеальной несжимаемой, лишенной капиллярных свойств, невесомой жидкости под действием постоянного перепада давления. Приводится численное решение задачи Релея. Под действием постоянного
перепада давления пузырь с исходного радиуса Ro=Rm замыкается полностью за время t* =0.91468 «іАрір. В этой задаче распределение
потенциала и скоростей на границе пузыря известно во все моменты времени. Если в начальный момент времени взять отрицательный перепад давления и соответствующее ему значение потенциала, то пузырь с
начального малого радиуса R0=QARm (/0 =0.00155 — J'Ap'/p) под
R»,
действием растягивающих напряжений увеличится до максимального
радиуса Rm) а затем при смене растягивающих напряжений избыточным
давлением перейдет в фазу замыкания. В этом случае полное время жизни
пузыря составит 2t*. Эта задача имеет существенно нелинейный характер и является прекрасной возможностью для тестирования численно алгоритма, отладки методики нахождения поля скоростей и подбора оптимальных параметров в алгоритме автоматического выбора шага по времени.
Также во второй главе исследуется задача об эволюции парогазового пузыря в безграничной идеальной несжимаемой невесомой, лишенной капиллярных свойств жидкости [54]. Данная задача отличается от задачи Релея предположением о том, что в начальный момент времени в пузыре находится некоторое количество газа, который в дальнейшем ведет себя адиабатически. Наличие газа в пузыре оказывает принципиальное влияние на его поведение. Как бы ни было мало начальное давление газа, при замыкании пузыря из начального радиуса R0 = Rm, давление в пузыре увеличивается, приостанавливая рано или поздно процесс замыкания пузыря. Пузырь замыкается до минимального радиуса, зависящего от параметра газосодержания (такая оценка приводится в работе [74]), после чего вновь расширяется до исходного радиуса, оставаясь сферическим. Такие пульсации пузыря, то есть сменяющие друг друга этапы расширения и замыкания, могут продолжаться бесконечно, так как потерь энергии не происходит, потенциальная и кинетическая энергия переходят одна в другую. В процессе же численного моделирования удается пройти всего несколько пульсаций пузыря, после чего происходит разрушение численно счета из-за нарушения сферической симметрии пузыря. В работе [101] описывается численное моделирование пульсаций пространственного парогазового пузыря в безграничной идеальной несжимаемой жидкости с учетом поверхностного натяжения в поле переменного давления. Моделирование
осуществляется методом обобщенных вихрей (generalized vortex method), в этом случае удается пройти всего лишь несколько пульсаций пузыря, и что интересно потеря сферической симметрии пузыря в работе [101] происходит аналогично тому, что наблюдается в диссертационном исследовании.
Во второй главе, кроме описанных задач, исследуется консервативность МГЭ путем контроля за выполнением закона сохранения энергии, приведенного в статье [25],
В третьей главе рассматривается эволюция пузыря в безграничной жидкости под действие силы тяжести, а также в невесомой жидкости около плоских и угловых твердых стенок. Влияние этих факторов приводит к нарушению сферической симметрии пузыря, что особенно заметно па этапе замыкания, когда происходит развитие струи жидкости, внедряющейся в пузырь. Особое внимание уделяется анализу характеристик струйки, методика выделения струйки и расчета ее характеристик описывается в первом параграфе данной главы. В качестве основных характеристик струйки используются высота струйки и скорость ее вершины, а так же, кинетическая энергия струйки и всего пузыря. Для оценки возможности урона, наносимого струйкой твердой стенке, на основании рассуждений приведенных в книге В.К.Кедринского [47] строится характеристика — коэффициент пробоя струйки. Этот коэффициент позволяет оценить возможность нанесения урона твердой стенке, при условии, что струйка развивает необходимую для этого скорость.
В третьей главе диссертации так же изучается влияние силы тяжести и газосодержания на процесс эволюции пузыря. В случае наличия силы тяжести струйка развивается в направлении противоположном ее действию, высота и скорость струйки зависят от величины коэффициента плавучести. Значения коэффициентов плавучести берутся из описанных экспериментальных работ [79, 80, 90]. В этом же параграфе рассматривается влияние содержания газа на динамику пузыря в случае, когда пузырь утрачивает сферическою форму.
При эволюции пузыря около твердой стенки наблюдается описанный, как в экспериментальных, так и в посвященных численному моделированию работах, эффект притяжения пузыря твердой стенкой и последующего формирования кумулятивной струйки развивающейся по направлению к ней. Этот эффект наблюдается независимо от наклона твердой стенки и определяется только расстоянием до нее. Вопрос о том, на каком расстоянии от стенки должен находиться пузырь, чтобы ее влияние оставалось
заметным, также исследуется в третьей главе. Значение параметра газосодержания берется из описанных экспериментальных работ [84], Приводится сравнение вычисленных скоростей струйки со скоростями, описанными в экспериментальных работах [80], а также со скоростями, полученными при расчете осесимметричной задачи, описанной в работе [82]. Здесь же, проводится сравнение форм пузыря, полученных при решении задачи в пространственной постановке с задачами в осесимметричной постановке, реализованой в работе К.Е.Афанасьева и А.М.Гудова [8].
Кроме того, в третьей главе диссертации рассматривается вопрос о влиянии поверхностного натяжения на процесс эволюции пузыря. Рассматривается методика вычисления поверхностного натяжения в узлах сетки, являющихся стыком плоских треугольных элементов. Считается, что действие поверхностного натяжения существенно только для паровых кавитационных пузырей, содержание газа в которых пренебрежимо мало [48, 93], поэтому в работе приводятся расчеты с коэффициентами поверхностного натяжения для кавитационных пузырей. Силы поверхностного натяжения стремятся сжать каверну, образовавшуюся в жидкости [48], а также препятствуют деформации пузыря на этапе его замыкания, заметно сокращая время существования пузыря. Наличие газа в пузыре препятствует его замыканию, отчасти компенсируя действие сил поверхностного натяжения, что подтверждается численными расчетами. Тем не менее, как показывают численные расчеты, под действие поверхностного натяжения даже пузырь со значительным содержанием газа замыкается, оставаясь значительно более гладким, чем в случае отсутствия сил поверхностного натяжения.
Четвертая глава диссертации посвящена исследованию динамики парогазового пузыря около наклонных твердых стенок, а также изучению возможности предсказания направления миграции пузыря и последующего образования струйки с использованием импульса Кельвина.
В ряде работ, посвященных эволюции парогазовых пузырей, описывается методика предсказания направления миграции пузыря и последующего образования струйки с использованием величины, которую можно получить из вектора количества движения, сделав ряд допущений [76, 78, 82, 83]. В первом параграфе четвертой главы диссертации приводится последовательность рассуждений для получения данной оценки и обсуждается достоверность предсказаний, полученных с ее помощью.
Пузырь развивается из малого сферически симметричного пузыря, расположенного на расстоянии одного или двух максимальных радиусов от
твердой стенки. В процессе роста пузырь, как правило, сохраняет близкую к
сферической форму. В случае, когда стенка расположена на расстоянии
одного максимального радиуса от пузыря, при достижении максимального
объема пузырь несколько "уплощается" со стороны твердой стенки. На этапе
замыкания пузыря рассматриваются несколько вариантов дальнейшего
развития течения. Пузырь, замыкаясь, может образовывать острые грани или
вершины, например, пузырь принимает форму капли, что проводит к
дальнейшем разрушению численного расчета и предсказать дальнейшее
поведение пузыря в этом случае не представляется возможным. Другой
вариант развития событий, когда просматривается тенденция к разделению
пузыря на два отдельных более мелких пузыря. Наиболее же вероятный
вариант развития событий, когда на этапе замыкания формируется струйка,
внедряющаяся в пузырь. Характеристики струйки определяются
преимущественно положением твердой стенки и величиной коэффициента плавучести. В большинстве случаев расчет удается продлить до момента касания твердой стенкой противоположной стенки пузыря. Существует несколько как экспериментальных так и численных работ, описывающих дальнейшее преобразование такого пузыря в тороидальный и дальнейшую его эволюцию [81, 84].
При численном моделировании задачи эволюции пузыря все характеристики динамики пузыря получаются в безразмерном виде. Поэтому интересно перейти к размерным величинам и получить реальные характеристики различных пузырей. В данной работе рассматриваются кавитационные пузыри радиусом 1.27 и 10 мм, описанные в книге Р.Кнэппа, Дж. Дейли и Ф.Хэммита [48], кроме того, рассматривается пузырь, образующийся при взрыве мины с начальным радиусом 1м, описанный в книге ГЛамба [57], и пузырь, образуемый ядерным взрывом "Вигвам" мощностью ЗОКт, описанным в работе Дж.У. Притчетт [64]. Для всех этих пузырей приводятся размерные высоты и скорости струек, а также полное время эволюции пузырей.
Таким образом на защиту выносятся следующие результаты:
Разработан численный алгоритм решения нестационарных пространственных задач динамики жидкости со свободными границами методом линейных граничных элементов. Отработан алгоритм автоматического выбора шага по времени.
Решены новые нестационарные задачи динамики пространственного парогазового пузыря:
о пульсации парогазового пузыря, минимальный радиус которого зависит от параметра газосодержания. Минимальные радиусы пузыря, полученные при численных расчетах, удовлетворяют существующим аналитическим оценкам,
о взаимодействие парогазового пузыря с наклонными и угловыми твердыми стенками.
Разработана методика выделения струйки и определения ее характеристик, таких как высота струйки, скорость ее вершины, кинетическая энергия. Введена величина, позволяющая оценивать возможность урона, наносимого струйкой мишени, — коэффициент пробоя струйки.
Изучено в пространственной постановке влияние силы тяжести, содержания газа, поверхностного натяжения и близости стенки на процесс эволюции пузыря.
Исследуется консервативность численного алгоритма на основе выполнения закона сохранения энергии. Приводятся формы пузыря при его эволюции около различных плоских наклонных и угловых стенок на различном расстоянии от них, для различных значений параметра плавучести.
Приводится анализ направления развития кумулятивной струйки и возможного урона, наносимого струйкой твердой стенке в зависимости
от положения стенки и значения параметра плавучести, для различных типов пузырей.
Результаты работы по мере их получения докладывались и обсуждались на VI Международной конференции «Вычислительные методы в задачах волновой гидродинамики» (Новосибирск, 1998г.), Сибирской школе-семинаре «Математические проблемы механики сплошных сред» (Новосибирск, 1998г.), Международной научной конференции «Моделирование, вычисления, проектирование в условиях неопределенности - 2000» (Уфа, 2000г.), VII Международной конференции «Вычислительные методы в задачах волновой гидродинамики» (Новосибирск, 2000г.), международной научной конференции, посвященной 90-летию со дня рождения Г.Г. Тумашева (Казань, 2000г), конференции «Конференция молодых ученых, посвященная 10-летию ИВТ СО РАН» (Новосибирск, 2000), конференции «Молодые ученые Кузбассу. Взгляд в XXI век» (Кемерово, 2001), международной летней научной школы «Гидродинамика больших скоростей» (Чебоксары, 2002), а также регулярно на научном семинаре «Численные методы решения задач механики сплошной среды»
Кемеровского государственного университета (Кемерово, 1997-2003гг.).
Основные результаты диссертации опубликованы в статьях [6, 11, 32,
33, 35, 36, 37, 74], и тезисах и аннотациях докладов [4, 5, 7, 31, 34].
Автор диссертации выражает глубокую признательность научному
руководителю К.Е.Афанасьеву - за постоянное внимание, терпение,
творческие идеи и неоценимую помощь при выполнении работы.
Граничные элементы и интерполирующие функции
Также во второй главе исследуется задача об эволюции парогазового пузыря в безграничной идеальной несжимаемой невесомой, лишенной капиллярных свойств жидкости [54]. Данная задача отличается от задачи Релея предположением о том, что в начальный момент времени в пузыре находится некоторое количество газа, который в дальнейшем ведет себя адиабатически. Наличие газа в пузыре оказывает принципиальное влияние на его поведение. Как бы ни было мало начальное давление газа, при замыкании пузыря из начального радиуса R0 = Rm, давление в пузыре увеличивается, приостанавливая рано или поздно процесс замыкания пузыря. Пузырь замыкается до минимального радиуса, зависящего от параметра газосодержания (такая оценка приводится в работе [74]), после чего вновь расширяется до исходного радиуса, оставаясь сферическим. Такие пульсации пузыря, то есть сменяющие друг друга этапы расширения и замыкания, могут продолжаться бесконечно, так как потерь энергии не происходит, потенциальная и кинетическая энергия переходят одна в другую. В процессе же численного моделирования удается пройти всего несколько пульсаций пузыря, после чего происходит разрушение численно счета из-за нарушения сферической симметрии пузыря. В работе [101] описывается численное моделирование пульсаций пространственного парогазового пузыря в безграничной идеальной несжимаемой жидкости с учетом поверхностного натяжения в поле переменного давления. Моделирование осуществляется методом обобщенных вихрей (generalized vortex method), в этом случае удается пройти всего лишь несколько пульсаций пузыря, и что интересно потеря сферической симметрии пузыря в работе [101] происходит аналогично тому, что наблюдается в диссертационном исследовании.
Во второй главе, кроме описанных задач, исследуется консервативность МГЭ путем контроля за выполнением закона сохранения энергии, приведенного в статье [25],
В третьей главе рассматривается эволюция пузыря в безграничной жидкости под действие силы тяжести, а также в невесомой жидкости около плоских и угловых твердых стенок. Влияние этих факторов приводит к нарушению сферической симметрии пузыря, что особенно заметно па этапе замыкания, когда происходит развитие струи жидкости, внедряющейся в пузырь. Особое внимание уделяется анализу характеристик струйки, методика выделения струйки и расчета ее характеристик описывается в первом параграфе данной главы. В качестве основных характеристик струйки используются высота струйки и скорость ее вершины, а так же, кинетическая энергия струйки и всего пузыря. Для оценки возможности урона, наносимого струйкой твердой стенке, на основании рассуждений приведенных в книге В.К.Кедринского [47] строится характеристика — коэффициент пробоя струйки. Этот коэффициент позволяет оценить возможность нанесения урона твердой стенке, при условии, что струйка развивает необходимую для этого скорость.
В третьей главе диссертации так же изучается влияние силы тяжести и газосодержания на процесс эволюции пузыря. В случае наличия силы тяжести струйка развивается в направлении противоположном ее действию, высота и скорость струйки зависят от величины коэффициента плавучести. Значения коэффициентов плавучести берутся из описанных экспериментальных работ [79, 80, 90]. В этом же параграфе рассматривается влияние содержания газа на динамику пузыря в случае, когда пузырь утрачивает сферическою форму.
При эволюции пузыря около твердой стенки наблюдается описанный, как в экспериментальных, так и в посвященных численному моделированию работах, эффект притяжения пузыря твердой стенкой и последующего формирования кумулятивной струйки развивающейся по направлению к ней. Этот эффект наблюдается независимо от наклона твердой стенки и определяется только расстоянием до нее. Вопрос о том, на каком расстоянии от стенки должен находиться пузырь, чтобы ее влияние оставалось заметным, также исследуется в третьей главе. Значение параметра газосодержания берется из описанных экспериментальных работ [84], Приводится сравнение вычисленных скоростей струйки со скоростями, описанными в экспериментальных работах [80], а также со скоростями, полученными при расчете осесимметричной задачи, описанной в работе [82]. Здесь же, проводится сравнение форм пузыря, полученных при решении задачи в пространственной постановке с задачами в осесимметричной постановке, реализованой в работе К.Е.Афанасьева и А.М.Гудова [8].
Кроме того, в третьей главе диссертации рассматривается вопрос о влиянии поверхностного натяжения на процесс эволюции пузыря. Рассматривается методика вычисления поверхностного натяжения в узлах сетки, являющихся стыком плоских треугольных элементов. Считается, что действие поверхностного натяжения существенно только для паровых кавитационных пузырей, содержание газа в которых пренебрежимо мало [48, 93], поэтому в работе приводятся расчеты с коэффициентами поверхностного натяжения для кавитационных пузырей. Силы поверхностного натяжения стремятся сжать каверну, образовавшуюся в жидкости [48], а также препятствуют деформации пузыря на этапе его замыкания, заметно сокращая время существования пузыря. Наличие газа в пузыре препятствует его замыканию, отчасти компенсируя действие сил поверхностного натяжения, что подтверждается численными расчетами. Тем не менее, как показывают численные расчеты, под действие поверхностного натяжения даже пузырь со значительным содержанием газа замыкается, оставаясь значительно более гладким, чем в случае отсутствия сил поверхностного натяжения.
Четвертая глава диссертации посвящена исследованию динамики парогазового пузыря около наклонных твердых стенок, а также изучению возможности предсказания направления миграции пузыря и последующего образования струйки с использованием импульса Кельвина.
В ряде работ, посвященных эволюции парогазовых пузырей, описывается методика предсказания направления миграции пузыря и последующего образования струйки с использованием величины, которую можно получить из вектора количества движения, сделав ряд допущений [76, 78, 82, 83]. В первом параграфе четвертой главы диссертации приводится последовательность рассуждений для получения данной оценки и обсуждается достоверность предсказаний, полученных с ее помощью.
Пузырь развивается из малого сферически симметричного пузыря, расположенного на расстоянии одного или двух максимальных радиусов от твердой стенки. В процессе роста пузырь, как правило, сохраняет близкую к сферической форму. В случае, когда стенка расположена на расстоянии одного максимального радиуса от пузыря, при достижении максимального объема пузырь несколько "уплощается" со стороны твердой стенки. На этапе замыкания пузыря рассматриваются несколько вариантов дальнейшего развития течения. Пузырь, замыкаясь, может образовывать острые грани или вершины, например, пузырь принимает форму капли, что проводит к дальнейшем разрушению численного расчета и предсказать дальнейшее поведение пузыря в этом случае не представляется возможным. Другой вариант развития событий, когда просматривается тенденция к разделению пузыря на два отдельных более мелких пузыря. Наиболее же вероятный вариант развития событий, когда на этапе замыкания формируется струйка, внедряющаяся в пузырь. Характеристики струйки определяются преимущественно положением твердой стенки и величиной коэффициента плавучести. В большинстве случаев расчет удается продлить до момента касания твердой стенкой противоположной стенки пузыря. Существует несколько как экспериментальных так и численных работ, описывающих дальнейшее преобразование такого пузыря в тороидальный и дальнейшую его эволюцию [81, 84].
Замыкание газонаполненного пузыря
Второй подход, как уже упоминалось, разработан специально для построения сетки на поверхности сферы (эллипсоида). Исходным приближением поверхности сферы является икосаэдр, состоящий из 12 узлов и 20 элементов, все узлы которого лежат на поверхности дискретизируемой сферы. На каждой итерации каждый исходный поверхностный треугольник разбивается на четыре следующим образом: каждая его сторона делится пополам новым узлом, который сдвигается в радиальном направлении на поверхность сферы, после чего все новые узлы соединяются ребрами. Количество элементов при переходе к новому уровню дискретизации увеличивается в четыре раза и может быть вычислено по формуле М - 5-І7" (п - номер уровня дискетизации, первый уровень - икосаэдр), а количество узлов N = 5-22" ] +4. Очевидно, что количество элементов сетки растет очень быстро, поэтому реально удается получить од ну-две приемлемые расчетных сетки. Количество узлов и элементов для каждого уровня дискретизации указано в таблице 1 , первые уровни дискретизации показаны на рис.5. Подобный, но немного усложненный подход, подход к построению поверхностной сетки описан в работе [111].
Преимущество второго метода состоит в том, что построенная с его помощью сетка является более равномерной, что подтверждается графиками, приведенными на рис,6, которые отражают зависимость отношения наибольшей площади элемента к наименыпей в зависимости от числа элементов для первого и второго типа сеток. Кроме того, сетка, полученная вторым способом, является более однородной, так как для такой сетки каждый узел окружен шестью элементами, за исключением 12 узлов, принадлежащих исходному икосаэдру, которые окружают по пять элементов. В сетках первого типа все узлы также окружены шестью элементами, но точки, являющиеся узловыми точками опорных зон, может окружать произвольное число элементов в зависимости от расположения опорных зон. В нашем случае такие точки окружены четырьмя или шестью элементами, но избежать появления узлов, окруженных четырьмя элементами, не удается.
Для решения полученной линейной задачи на каждом временном шаге был выбран метод граничных элементов. Его суть состоит в преобразовании дифференциального уравнения в частных производных, которое описывает поведение неизвестных функций внутри и на границе области, в интегральное уравнение, записанное по границе [14], Перейти к интегральной формулировке в нашем случае можно используя третью формулу Грина где р - функция потенциала поля скоростей, q = dq jdn - его нормальная производная, п - внешняя по отношению к области нормаль (в данном случае направленная внутрь пузыря), р - фундаментальное решение уравнения Лапласа, которое в пространственном случае записывается в виде р" = 1. 4лг(, ;), г( , ) - растояние между точками х и %, % -текущая точка на границе области Г, х будем называть точкой наблюдения, х также находится на границе Г. Коэффициент С = со 2тг, со представляет собой телесный угол, под которым видна поверхность из точки X. Реализация метода граничных элементов состоит из следующих этапов: 1. Граница области, в данном случае Г, разбивается на ряд элементов, на которых предполагается, что потенциал и его нормальная производная изменяются в соответствии с выбранными интерполирующими функциями. На каждом элементе в соответствии с порядком аппроксимации задается необходимое число узлов. 2. Для каждой узловой точки, которая в данном случае будет называться точкой наблюдения, записывается дискретная форма уравнения (1.18), связывающего потенциал и его нормальную производную в каждом узле. 3. Интегралы по каждому элементу вычисляются аналитически или с помощью схем численного интегрирования. 4. После наложения заданных граничных условий получается система линейных алгебраических уравнений, относительно остальных неизвестных значений потенциала или его нормальной производной в узлах, которая решается прямым или итерационным методом. Решение системы дает эти неизвестные значения потенциала и его нормальной производной в узлах границы. 5. По известным значениям потеїщиала и его нормальной производной в узловых точках находится распределение поля скоростей на границе области, что дает возможность перейти к следующему шагу по времени. В случае необходимости значения функции потенциала внутри области могут быть найдены по известным значениям потенциала и его нормальной производной на границе с помощью численного интегрирования выражения (1.18).
Сравнительный анализ расчета с использованием пространственной и осесимметричной моделей
Очевидно, что результирующая система алгебраических уравнений будет зависеть от свойств интегральных уравнений, записанных для каждого узла на границе области.
В нашем случае имеем задачу Дирихле и исходный дифференциальный оператор сводится к интегральному уравнению Фредгольма I рода вида которое является некорректно поставленным по Адамару [70]. Задача считается корректной по Адамару, если выполняются условия существования решения, его единственности и условие, что малым ошибкам исходных данных соответствуют малые ошибки решения (условие устойчивости решения). Задача решения уравнения (1.36) является некорректной, в первую очередь, из-за нарушения условия устойчивости решения, когда даже очень малые относительные ошибки правой части f{x), а также ядра K{x,s) И метода решения могут
приводить к таким ошибкам, что численное решение не будет иметь практически ничего общего с точным [23]. В нашем случае дискретизация границы, аппроксимация потенциала и его нормальной производной линейными функциями, а также численное вычисление интегралов вносят в интегральное уравнение определенную погрешность. Кроме того, при вычислении диагональных коэффициентов матрицы G в случае, когда точка наблюдения является одной из вершин элементов, по которым происходит интегрирование, приходится иметь дело с интегралами, существующими только в смысле главного значения по Коши. Решение задачи (1.36) может терять устойчивость вследствии плохой обусловленности результирующей СЛАУ при достаточно большом числе элементов. Поэтому, для численного решения таких уравнений должны использоваться методы регуляризации некорректных задач. В качестве одного из таких подходов можно привести регуляризующий алгоритм Тихонова [69], который основан на переходе к вспомогательному уравнению Фредгольма II рода (уравнению Эйлера для сглаживающего функционала). Такой алгоритм обеспечивает устойчивость численного решения к погрешностям входной информации, но является алгоритмически значительно более сложным при программной реализации. В нашем случае в использовании такого алгоритма нет необходимости, так как квадратурные формулы для решения уравнений Фредгольма I рода с ядрами типа функций Грина порождают саморегуляризирующие алгоритмы, в которых параметром регуляризации является шаг квадратурной формулы [29]. Если же интегральные коэффициенты имеют сильную особенность, как в нашем случае, которая вычисляется в смысле главного значения по Коши, то при решении результирующей системы линейных алгебраических уравнений также наблюдается эффект "саморегуляризации" [58].
На каждом шаге по времени свойства матрицы неизвестны, а их выяснение требует значительных временных затрат. Только факт самосопряженности установить легко, свойства знакоопределенности и неособенности установить достаточно сложно. Очень полезной характеристикой в этом смысле является число обусловленности, которое выполняет роль коэффициента увеличения относительной ошибки [46]. Другими словами число обусловленности cond(A) является мерой относительного расстояния от Л до множества вырожденных матриц. В задаче решения системы линейных уравнений если cond{A)&\q, — машинная точность, то метод исключения Гаусса дает решение с d-q правильными машинными знаками [71, 72]. В нашем случае число обусловленности не достигает таких больших величин как можно было бы ожидать (не более 105), что при большом числе доступных машинных разрядов (не меньше 16) позволяет получить приемлемую точность решения. Показательно, что увеличение машинной точности до 10 32, что возможно при использовании Фортрана, поставляемого с современными версиями операционной системы Solaris, не позволяет добиться повышения точности счета. При расчете тестовых задач были использованы как точные методы решения СЛАУ (метод Гаусса с выбором главного элемента), так и итерационные (метод Гаусса-Зейделя, метод нелинейной регуляризации (В.Н.Трушникова)) [10, 12], наиболее же приемлемым оказался метод Гаусса с выбором главного элемента с последующим итерационным уточнением из библиотеки IMSL Microsoft Fortran Power Station. В данном случае точность решения методом Гаусса удается повысить , выполнив его итерационное уточнение, которое при сравнительно небольшом числе обусловленности матрицы А позволяет выполнить вычисления с относительной точночтыо сравнимой с єм. Его идея заключается в следующем. Пусть х w & - соответственно точное и приближенное решения системы линейных уравнений. Пусть х = +5, где d-(S1,S2,...,Sn)T - поправка решения J, a r = {rur2,...,rn)T - невязка для приближенного решения : г = Ь-А. Тогда имеем A(x±S)=b или AS = г. То есть, чтобы найти уточнение 5, следует решить систему A3-г с прежней матрицей А. Для ее решения используется раннее найденное разложение матрицы А на верхне- и нижнетреугольную матрицы, чем и обусловлены сравнительно небольшие временные затраты на итерационное уточнение [15].
Характеристики струйки. Переход к размерным величинам
Очевидно, что схлопывание пузыря при описанных условиях -существенно нелинейный процесс. При г»0 скорость стремиться к бесконечности как Г3,2. Известное аналитическое решение хорошо согласуется с результатами экспериментов [48], за исключением последних моментов существования пузыря, когда скорость становится бесконечной. Аналитическое решение дает возможность выбрать оптимальные параметры движения по времени Дгтіп, Агтах и %, верифицировать работу численного алгоритма в целом, и в особенности, методики нахождение скоростей и движения по времени. В результате тестирования максимальный и минимальный шаг по времени были определены как Armax=0,01, Дгтіп =0.0001, =0.02. Для численного решения параметры имеют значения а = 0.0, /? = 0.0, = 0.0, в начальный момент времени пузырь представляет собой сферу SQ радиуса l.0tfmax с центром в начале координат, начальное значение потенциала )(0); =0. Рисунок 15, на котором приведено сравнение решений с начального момента времени до момента коллапса пузыря 0.0 г г показывает хорошее соответствие решений. Рисунок 16 демонстрирует точность вычисления значения нормальной скорости q в моменты близкие к коллапсу пузыря 0.83 г 0.92, для аналитического решения момент коллапса пузыря и момент приходящийся на шаг по времени перед ним опущены, так как значение нормальной скорости растет чрезвычайно быстро и к моменту г обращается в бесконечность. Рисунок 16 демонстрирует сравнение численного и аналитического значений нормальной скорости не с начального момента времени, а с г = 0.83, так как графики, показанные с начального момента, практически совпадают. Согласно аналитическому решению пузырь должен сохранять сферическую симметрию все время существования, форма пузыря в процессе его эволюции показана в таблице 3. Теперь предположим, что процессу замыкания пузыря под действием избыточных давлений предшествует его рост под действием растягивающих напряжений. В этом случае давление паров в пузыре должно быть больше гидростатического давления Aps = рн—р00 0. Тогда уравнение (2.1) примет вид А интегрирование уравнения (2.4) с использованием начальных условий г=0 =0 дает следующее выражение для скорости Пузырь расширяется с нулевого радиуса и значение радиуса как бы движется в обратном направлении по графику радиуса пузыря от времени в случае его замыкания. По достижении максимального радиуса =1 растягивающие напряжения Aps сменяются избыточными давлениями Ар и пузырь переходит в фазу замыкания, описаную ранее. Полное время эволюции пузыря в этом случае составит = 2г = 1.82936. Так как начать расчет с нулевого радиуса не представляется возможным, в качестве начальных условий возьмем пузырь малого радиуса R0=0ARm (г =0.00155), в качестве начального возьмем значение потенциала из уравнения (2.4). В этом случае зависимость радиуса от времени перестает быть монотонной. Более того, в следствии «склейки» решений график зависимости радиуса от времени образует «излом» при т = т", когда растягивающие напряжения сменяются сжимающим давлением. На рисунке 17 показаны аналитическое и численные решения, полученные при подборе оптимальных значений Атт[п, Дгтах и х- Также на этом рисунке приводятся увеличенные моменты достижения пузырем максимального радиуса и последние моменты перед замыканием. Из рисунка видно, что в случае численного моделирования процесса расширения и последующего замыкания пузыря график зависимости радиуса пузыря от времени получается гладким, В остальном численное решение показывает хорошее соответствие аналитическому, сходимость алгоритма в процессе подбора шага очевидна. Для подбора шага было перебрано значительно большее число наборов параметров Armin, Дгтах и Z, но для наглядности рисунка приводится всего пять графиков численных решений для различных сочетаний этих параметров.
Принципиальное влияние на характеристики движения пузыря оказывает наличие в нем газа. Начальное давление газа, пусть даже небольшое, при уменьшении объема пузыря повышается, противодействуя движущейся границе, таким образом, процесс замыкания рано или поздно будет приостановлен и сменится процессом расширения. Исследуем характеристики влияния газа, для чего в р„ включим слагаемое Po{VQ/V)y, где V - объем пузыря и V(0) = V0 - объем пузыря в начальный момент времени, предполагая, что сжатие газа происходит по адиабатическому закону. Тогда уравнение движения границы пузыря примет вид [54] у - отношение удельных теплоємкостей газа. Приняв для упрощения у=А!Ъ вычислим первый интеграл уравнения (2.4) [54] Таким образом, пузырь, а начальный момент времени представляющий собой сферу 0 радиуса 1Лет, сохраняя сферическую форму сжимается до минимального радиуса #min(/?), после чего вновь расширяется до исходного радиуса. Так как данная модель не предусматривает потерь энергии, то пульсации пузыря могут продолжаться бесконечно. При численном моделировании удается пройти всего несколько пульсаций, после чего происходит разрушение процесса численного счета из-за нарушения сферической симметрии пузыря. В таблице 4 для значений р от 0.4 до 0.9 приведены аналитическое Rminau. и численные значения минимальных радиусов для двух первых пульсаций пузыря, Я1ттчисл. и Я ачисл., соответственно, полное время жизни пузыря до разрушения расчета теУ полная энергия Е и колебание полной энергии АЕ в процентах. Чем больше р, тем больше частота и меньше амплитуда пульсаций. Очевидно, что из-за малой амплитуды наибольшее количество пульсаций удается пройти при р=0.9 — почти пять, в то время как при P=0J - три неполных пульсации, а при р=0А всего две неполных пульсации. Интересно, что за исключением расчета Р 0.6 разрушение расчета происходит на этапе роста пузыря. На рис 17 для всех значении р приведены по три формы пузыря, когда он последний раз достигает максимальных и минимальных объемов, а так же последний момент перед разрушением расчета.
