CC BY
97
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА
Физико-математические пауки
УДК 519.6:532:533
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОНТАКТНОГО ВЗАИМОДЕЙСТВИЯ СЖИМАЕМЫХ СРЕД НА ЭЙЛЕРОВЫХ СЕТКАХ
А. А. Агапип, Т. С. Гусева
Аннотация
Представлен краткий обзор работ по численному моделированию задач контактного взаимодействия сжимаемых сред с большой разницей акустических импедапсов при сильных деформациях контактной границы и при наличии в средах ударных воли. Приведены основные положения реализованной авторами методики расчета таких задач па основе метода CIP-CUP (Constrained Interpolation Profile Combined Unified Procedure), в котором применяются эйлеровы сетки без явного выделения контактных границ, а исходные уравнения движения сжимаемой жидкости записываются в терминах пекопсервативпых переменных. Представлены результаты расчетов тестовых задач, имеющих аналитическое решение, подтверждающие работоспособность созданного алгоритма. Возможности методики проиллюстрированы па задачах об ударе осесимметричпой высокоскоростной струи жидкости по жесткой степке и по топкому слою жидкости па жесткой степке. Результаты расчетов удовлетворительно согласуются с известными численными решениями. полученными методом адаптивпо-подвижпых сеток с явным выделением межфазпой границы.
Ключевые слова: течения многофазных сред, большие деформации межфазпой границы. эйлерова сетка, «улавливание» контактной границы, метод CIP-CUP.
Введение
Многие важные с практической точки зрения явления характеризуются одновременным наличием сильных ударных волн и контактных границ, которые претерпевают большие перемещения, деформации, нарушения топологии. Примеры таких явлений подводные взрывы вблизи свободной поверхности, взаимодействие струй, капель, пузырьков с ударными волнами, твердыми стенками, границами раздела жидкостей. Для описания подобных явлений необходимо, как правило, применять уравнения динамики сплошных сред с учетом эффектов вязкости, теплопроводности н дополнять их сложными уравнениями состояния. Решение задач при таком подходе приходится находить численно. Обзоры большого количества используемых в настоящее время методов численного решения задач динамики сплошных сред при наличии в них сильных ударных волн, разнообразных контактных границ можно найти в [1 11]. Имеющиеся численные методы принято условно разбивать на лагранжевы. произвольные эйлерово-лагранжевы и эйлеровы.
В лагранжевых методах ячейки расчетной сетки являются частицами среды. При этом контактные границы всегда совпадают с координатными поверхностями и естественным образом деформируются в соответствии с деформацией частиц, так что никаких специальных процедур для описания движения контактных границ не требуется. Поэтому для задач с изменяющимися контактными границами лагранжевы методы являются наиболее эффективными. Апироксимационные соотношения лагранжевых методов относительно просты, поскольку в лаграижевой форме уравнений динамики сплошных сред отсутствуют конвективные слагаемые.
Однако в силу недивергентного характера лагранжевой формы уравнений динамики сплошных сред лагранжевы методы неконсервативны. Это означает, что при исследовании задач с интенсивными ударными волнами нужно вводить искусственную вязкость [12. 13]. Кроме того, при сильных деформациях частиц среды, например. в потоках с формирующейся завихренностью, лагранжева сетка сильно искажается, что приводит к потере точности и росту погрешностей. В таких случаях приходится проводить редискретизацию изучаемых сред, другими словами, перестройку расчетной сетки, и выполнять интерполяцию параметров ячеек со старой сетки на новую. При этом преимущества лагранжевых методов сохраняются. но только до тех пор. пока процедуры перестройки сеток и интерполяции относительно просты и выполняются нечасто.
Проблемы потери точности расчетов, обусловленные деформацией вычислительных сеток, значительно уменьшаются при использовании произвольных лагранжево-эйлеровых (ALE. Arbitrary Lagrarigiari-Eulerian) методов (или методов подвижных сеток). В этих методах существенные границы задачи (ударные волны, контактные границы, свободные поверхности, жесткие стенки и т. д.) всегда совпадают с координатными поверхностями (или их частями). При этом нормальные составляющие скорости перемещения узлов расчетной сетки, расположенных на границе, совпадают со скоростью самой границы, а в касательном направлении эти узлы могут перемещаться произвольным образом. Перемещение остальных узлов сетки, не совпадающих с границей, может быть произвольным во всех направлениях. независимо от движения частиц среды. Эту свободу можно использовать для повышения экономичности расчетов. Так. можно строить адаптивные сетки со сгущением в областях с большими градиентами решения и разрежением в областях с гладким поведением решения, сетки с плавным изменением характерных размеров соседних ячеек, сетки с сохранением четырехугольных ячеек, близких к квадратным или прямоугольным и т. д. При этом необходимо вводить в математическую постановку задачи дополнительные уравнения движения узлов сетки. Исходные уравнения в произвольных лагранжево-эйлеровых методах можно представить как в недивергентной, так и в дивергентной формах. Поэтому расчет интенсивных ударных волн в рамках этих методов можно проводить без введения искусственной вязкости, воспользовавшись для этого консервативными схемами. Если лагранжевы методы считаются оптимальными в задачах с относительно малыми деформациями самих сред, то произвольные лагранжево-эйлеровы методы оптимальны при относительно небольших деформациях существенных границ задачи. под которые подстраиваются используемые в этих методах расчетные сетки. Однако если существенные границы сильно деформируются, возникают, исчезают или меняют свою связность, то при использовании произвольных лагранжево-эйровых методов также приходится выполнять перестройку сетки и проводить соответствующую интерполяцию численного решения на новую сетку (иногда с выполнением законов сохранения массы, импульса, полной энергии). Естественно, что все это усложняет алгоритмы расчета и снижает их экономичность пропорционально частоте применения процедур перестройки сетки и интерполяции решения.
В эйлеровых методах используются стационарные (декартовы или криволинейные) сетки. При их использовании не возникает проблем, связанных с перестройкой сеток и интерполяцией численного решения с сетки на сетку. Как и схемы произвольных лагранжево-эйлеровых методов, эйлеровы схемы могут быть как консервативными. так и неконсервативными, что зависит от исходной математической формулировки задачи (она может быть как дивергентной, так и недивергентной) и способа ее аппроксимации. Для расчета задач с сильными ударными волнами широко применяются консервативные методы без явного выделения газодинамических разрывов, в которых не нужно использовать искусственную вязкость.
Наиболее популярными среди них является метод Годунова [14] и его модификации повышенного порядка точности (повышение точности достигается за счет применения эффективных ограничителей производных).
Существующие методы описания перемещения контактной границы по эйлеровой сетке можно условно подразделить на две большие группы. Первую группу составляют методы, в которых подвижные контактные границы выделяются явно как совокупность поверхностных ячеек или узлов [10. 15 20]. В англоязычной литературе они называются FT-мстодами (front-tracking methods), а в терминологии работы [4] методами дискретных лагранжевых маркеров. Перемещение контактной границы (изменение позиции идентифицирующих ее маркеров) осуществляется явно путем решения уравнения
dx
It =U'
где u - скорость элемента контактной границы, х - его радиус-вектор, t - время. Поверхностная лагранжева сетка маркеров контактной границы, размерность которой на единицу меньше пространственной размерности рассматриваемой задачи, движется по эйлеровой сетке. Определять параметры сред на межфазной границе (жидкость газ) в случае ее явного выделения можно, например, решая задачу Ри-мана о распаде произвольного разрыва с единым уравнением состояния для обеих сред двучленным уравнением [14] или универсальным уравнением "stiffened gas" [21]. В целом, FT-методы характеризуются относительно высокой эффективностью расчета движения контактных границ, но при изменении топологии границ нлн при обобщении на многомерный случай алгоритмы значительно усложняются.
Во второй группе методов расчета контактных границ на эйлеровой сетке подвижные контактные границы улавливаются неявно посредством пространственной функции у(х, t) (х - радиус-вектор точки пространства), идентифицирующей контактирующие среды. В англоязычной литературе такие методы называются FC-методами (front-capturing methods), а в терминологии работы [4] методами непрерывных лагранжевых маркеров. Функция-идентификатор у рассчитывается наряду с остальными характеристиками сред. Основная идея применения функции-идентификатора заключается в том, что ее значения сохраняются вдоль лагранжевых траекторий частиц среды. Изменение функции у подчиняется уравнению переноса
t+u-уу=0'
где и локальная скорость среды.
Существующие FC-методы различаются способом задания и физическим смыслом функции-идентификатора. Описание ряда широко распространенных вариантов FC-методов можно найти в [3, 22]. К ним, в частности, относятся методы Volume Of Fluid (VOF) и Level Set (LS). В методе VOF функция-идентификатор у имеет смысл объемной доли данной среды в ячейке фракционного объема. При этом каждая из контактирующих сред характеризуется своей функцией-идентификатором, которая имеет значение 1 в ячейках, полностью занятых со-
0
01
предложен для несжимаемых и слабосжнмаемых жидкостей [23]. Сейчас существуют модификации этого метода для задач с полным учетом сжимаемости [24, 25]. В рамках этого метода для повышения точности определения положения контактной границы внутри ячейки либо применяются специальные алгоритмы реконструкции границы на каждом временном шаге на основе рассчитанной функции-у
у
метода VOF усложняются при изменении топологии границы и при обобщении на многомерный случай. В методе LS функция-идентификатор имеет смысл знакопеременного расстояния до контактной границы (знак функции идентифицирует среды по разные стороны от границы контакта) [30. 31]. При численном решении уравнения переноса функция у может терять свойства расстояния, поэтому в рамках метода LS приходится применять специальные процедуры ее ренннцналнзацнн.
Методы расчета течений с контактными границами подразделяют также на два класса в зависимости от того, как трактуется в них контактная граница. К первому-классу относят методы, в которых контактная граница рассматривается как четкая поверхность (линия), по разные стороны которой находятся разные среды. В англоязычной литературе такие методы называются sharp interface methods (SIM). Второй класс составляют методы, в которых контактная граница представляется не как четкая поверхность, а как тонкий переходный слой, в котором характеристики среды (некой смеси контактирующих сред) плавно переходят от параметров одной среды к параметрам другой. В англоязычной литературе такие методы называются diffusive interface methods (DIM). К методам класса SOI относятся все FT-методы и ряд FC-методов. Так. усилия, прилагаемые в FC-методах VOF и LS для того, чтобы функция-идентификатор у сохраняла заявленный физический смысл (объемной доли жидкости в ячейке для VOF и знакопеременного расстояния до контактной границы для LS), нужны для установления как можно более точного положения четкой границы между несмешнвающнмнся средами. В силу этого методы VOF и LS относятся классу SIM. В FC-методах класса DIM, как и в методе VOF, функция-идентификатор у (называемая также "indicating function", "color function", "donsity function" или, по аналогии с VOF, "fractional volume") в областях, занятых разными средами, принимает различающиеся постоянные значения. В начальный момент времени она имеет ступенчатый профиль па границе контакта. В ходе расчетов вследствие численной диффузии этот профиль постепенно размывается. При этом точное положение контактной границы становится неопределенным, как и в технологиях расчета ударных волн без их явного выделения. В случае необходимости граница контакта идентифицируется либо как область больших градиентов функции-идентификатора, либо как изоиоверхность ее среднего значения. Таким образом, хотя контактная граница между несмешива-ющимися жидкостями является в реальности четкой, в рамках методов класса DIM она заменяется на тонкую переходную область, занятую смесыо, сформированной контактирующими средами. С измельчением сетки переходная область превращается в разрыв. Для переходной области должны быть определены соответствующие переходные значения параметров уравнения состояния или само переходное уравнение состояния для "смеси". При этом функция-идентификатор может интерпретироваться как эффективный показатель адиабаты [32, 33], некая функция от показателя адиабаты [34], эффективная постоянная двучленного уравнения состояния (например, модели "stiffonod gas") [7, 8], массовая [35 38] или объемная [3, 39 45] доля среды в "смеси".
Методы DIM с размыванием контактной границы являются с точки зрения описания ее положения наименее точными, но в то же время они гораздо более просты в применении. Они легко обобщаются на многомерный случай и без дополнительных сложностей позволяют определять положение контактирующих сред в случае переменной топологии занимаемых ими подобластей. Достоинством этих методов является еще и то, что обусловленная численной диффузией ширина переходной области в них может быть ограничена при использовании определенных методик [11, 43, 46 48].
В последнее время появилось много гибридных методов (см. ссылки в [9]), в той или иной мере объединяющих достоинства сочетаемых в них подходов: смешанные алгоритмы LS-VOF, particle-LS, marker-VOF, LS-FT и т. д.
Таким образом, основное преимущество DIM-методов с размыванием контактных границ перед методами SIM с явным выделением контактных границ это относительная простота алгоритмов описания границ сложной топологии, включая фрагментацию, коалесценцию. возникновение и исчезновение поверхностей. В силу этого DIM-методы выглядят весьма привлекательными. Однако их использование для расчета контактного взаимодействия сжимаемых сред в сочетании с консервативными эйлеровыми схемами, основанными на сквозном расчете разрывов. оказалось проблематичным [32]. Дело в том. что консервативные эйлеровы методы сквозного счета всегда сопровождаются численной диффузией разрывов на несколько расчетных узлов. Эта диффузия, весьма желательная при расчете ударных волн в однородных средах, вызывает проблемы в случае расчета течений с контактными границами. Так. численная диффузия приводит к размыванию профиля плотности на контактной границе. В результате возникает переходная зона с неопределенной средой, давление которой уже нельзя выражать через консервативные переменные (массу, импульс и полную энергию), применяя для этого одно из уравнений состояния контактирующих сред или какое-либо уравнение состояния "смеси", поскольку это приводит к возникновению нефизичных осцилляций давления. Эти осцилляции не связаны с порядком точности схемы, поскольку наблюдаются даже при использовании схем первого порядка точности [22]. Нужно отметить, что FT-методы не имеют такого недостатка, поскольку контактная граница в них выделяется явно без какого-либо размывания.
Избежать нефизичных осцилляций давления в рамках DIM-иодхода можно, воспользовавшись иекоисервативиыми схемами сквозного счета, основанными на применении уравнения эволюции давления [32. 49. 50]. Можно использовать и консервативные схемы, но с локальной неконсерватнвной коррекцией полной энергии в окрестности контактной границы также на основе уравнения эволюции давления [3, 22. 51].
С другой стороны, в методах класса SIM с четким определением контактной границы применение консервативных эйлеровых схем сквозного счета, например схемы Годунова, также сопровождается некоторыми трудностями. При движении контактной границы по эйлеровой сетке возникают дробные ячейки (cut cells). Они имеют нерегулярную форму, так что для их обработки требуются специальные процедуры. точность которых обычно ниже точности алгоритмов расчета регулярных ячеек [1. 52]. Кроме того, размеры дробных ячеек могут быть очень малыми, что приводит к жестким ограничениям на шаг по времени.
Избежать проблем с дробными ячейками в рамках SIM-иодхода позволяет идеология фиктивной жидкости GFM (Ghost Fluid Method) [53]. Согласно этой идеологии поверх ячеек с реальной средой на "другой стороне" контактной границы вводится слой фиктивных ячеек, заполненных фиктивной жидкостью, термодинамически подобной жидкости "с этой стороны". После независимого расчета каждой из ''однородных" сред (реальная среда с одной стороны контактной границы термодинамически эквивалентная ей фиктивная среда по другую сторону границы) в зависимости от положения границы контакта в узлах сохраняется только решение. соответствующее реальной среде. Для доопределения фиктивных ячеек, входящих в шаблон схемы расчета ячейки с контактной границей, и для модификации потоков на гранях такой ячейки используется решение задачи Римана о распаде разрыва на самой контактной границе. В случае контакта сред с большой разницей акустических импедансов. то есть контакта сред типа газ жидкость с сильно различающимися уравнениями состояния, метод фиктивной жидкости GFM потребовал ряд модификаций [54 56].
В классе SIM-методов можно отметить метод I-SIM [57]. в котором движение контактной границы и динамика контактирующих сред являются полностью
(нелинейно) связанными, что исключает возникновение погрешностей, обусловленных операторным расщеплением. Кроме того, метод I-SIM является универсальным в том смысле, что он пригоден для расчета как сжимаемых, так и несжимаемых течений.
Прорыв в разработке подобных методов, которые можно применять для расчета контактного взаимодействия сред при произвольных скоростях их движения с учетом и без учета сжимаемости, был связан с появлением метода CIP-CUP (Constraint Interpolation Profile Combined and Unified Procedure) [49]. Этот метод основан на расщеплении уравнений Эйлера в переменных плотность, скорость, давление на конвективную и неконвективную части. Для расчета конвективной части применяется метод CIP один из вариантов иолулагранжевых методов, характеризующийся монотонностью и малой диффузией. Для расчета неконвективной части используется идеология методов расчета течений несжимаемых сред на основе уравнения для давления, обобщенная на случай сжимаемости. Данная процедура может применяться для расчета течений с широким диапазоном проявления сжимаемости н с контактными границами со значительным перепадом акустического импеданса. Несмотря на то что метод CIP-CUP не является консервативным, в сочетании с искусственной вязкостью [58] он демонстрирует хорошие возможности при расчете многофазных течений с ударными волнами. Поэтому сам метод CIP-CUP и его модификации широко применяются многими исследователями [43 45. 59 64]. Метод, идеологически эквивалентный методу [49], используется в работе [5] в качестве предобуславлнвателя для консервативной схемы с целыо обеспечения корректности решения на контактной границе.
В настоящей работе реализована методика расчета задач контактного взаимодействия сжимаемых сред, основанная на DIM-подходе с неявным улавливанием межфазных границ на эйлеровой сетке и уравнениях динамики сред в неконсервативной форме с применением для их решения метода CIP-CUP [43]. Такой подход сочетает в себе относительную алгоритмическую простоту, универсальность и эффективность при расчете взаимодействия сжимаемых и несжимаемых сред с сильными деформациями контактных границ.
Рассматривается контактное взаимодействие жидкости и газа сред, сильно различающихся между собой термодинамически и по проявлению эффекта сжимаемости. При этом допускаются конечные деформации контактной границы и наличие сильных ударных волн в жидкости и газе. Примерами таких задач являются охлопывание газовой полости в жидкости при прохождении в ней интенсивной ударной волны, сильное несферическое сжатие кавитационного пузырька вблизи тела с образованием на поверхности пузырька кумулятивной струйки, высокоскоростной удар движущейся в газе капли жидкости о жесткую стенку и т. д.
Уравнения динамики взаимодействующих сред без учета эффектов вязкости и теплопроводности записываются относительно плотности д, скорости и и давления р в виде
1. Постановка задачи и основные уравнения
в V-u,
Vp
в
^p + u •Vp = -eCS V • u,
где Св - скорость звука, определяемая общим для обеих сред выражением
Св = у 7 (Р + в).
Для жидкости входящие в это выражение постоянные 7 и В являются константами уравнения состояния Тэта, а для газа В = 0, 7 - показатель адиабаты. В расчетах скорость звука вычисляется по формуле
Св = уСв,; + (1 - у)Св,д,
где С^,; - скорость звука в жидкости, Свд - скорость звука в газе, у - функция-идентификатор среды. В начале расчета она определяется следующим образом:
1
У(х) = \
10, в области, запятой газом,
где х - радиус-вектор точки пространства. Изменение функции у описывается уравнением переноса
^ + и • У у =0, (2)
где и локальная скорость среды. При численном решении этого уравнения функция-идентификатор может принимать промежуточные значения между 0 и 1.
В таком случае контактная граница находится в области ненулевого градиента у
Как и в работе [43], решение системы (1) находится расщеплением уравнений на конвективную и ноконвоктивную части. Конвективная часть аппроксимируется следующим образом:
в* - в"
At и* -и"
At
p* - p"
+ и" • Ve" = 0,
- + и" • Vu" = 0, (3)
+ и" • Vp" = 0,
At
Звездочкой обозначаются величины после расчета конвективной части. Неконвективная часть аппроксимируется соотношениями
в"+1 _ в*
в д в = -e*V^u "+1, At
Vp"+1 (4)
At в*
p А p = -в*С82^и "+1.
At s
Для расчета конвективной части применяется метод CIP (Constrained Interpolation Profile) [43] один из вариантов полулагранжового метода. В этом методе на каждом временном шаге требуется определить положение и характеристики лагранжовой частицы, которая придет в рассматриваемый узел эйлеровой сетки. В частности, решение уравнения
f + u •Vf = 0 (5)
Рис. 1. Сплошная кривая - интерполяционный профиль функции в момент Ь — ДЬ; штриховая - этот же профиль, смещенный в направлении движения среды на расстояние иДЬ
записывается в виде лагранжева инварианта f (х, t) = f (х',t — At) (рис. 1), где х' - координаты частицы среды в момент t — At, которая, перемещаясь со скоро-
t
'
'
f( ', t — At) f( ', t — At) '
строится интерполяционный профиль функции. Его построение является ключевым моментом в CIP-методах. В настоящей работе применяется вариант RCIP (Rational-Cubic Interpolation Propagation) [43], в котором интерполяционная функция в одномерном случае в г-й ячейке представляется в виде
Д1* (х) = I Е арврХр I 1х. (6)
\0<р<1 У 0<1х<3
Здесь X = х — хг,
«о = 1, во = 1, «1 = ^и дг дЫр < 0, «1 = 0 при дг дыр > 0, в1 = [ — дг)/(дгир — $)| — 1 ]Дг\
Со = /г, С = дг + /г«1в1, С = 5г«1в1 + ($г — дг)Дгг1 — СзДг,
Сз = [дг — $ + (д гир — $г)(1 + а^Дг)] Д-2,
где гир = г — 1 щи и > 0, гир = г + ^и и < 0, Дг = хгир — хг, дг = дх/г, $ =
= (/гир — /г)/Дг • Все эти коэффициенты вычисляются по известным на п-м вре-
/
д = дх/. Определенная таким образом интерполяционная функция Д10 (х) обладает свойством снижения порядка интерполяции в областях с резким изменением
а1 = 0
областях, где решение гладкое (то есть в1 ^ 0). В ячейках, где решение имеет локальный экстремум и а1 = 1, Д10(х) реализует интерполяцию более низкого порядка с сохранением выпуклости или вогнутости профиля решения. При этом решение уравнения переноса (5) имеет вид
/* = Д1° (хг — иД£) = С + ^ + ^ + ^
1 + «1в1£
где £ = —и. При вычислении £ в качестве и используется осредпеппое значение (иг + и**)/2, где и** - скорость в точке хг — игД£. Значение и** предварительно
рассчитывается также с помощью Т1С1Р при и = иг. В рамках Т1С1Р пропзвод-д/
1
не через узловые значения функции f, а из уравнения, полученного дифференцированием уравнения для f. Например, если f удовлетворяет уравнению
df + df G dt + u дх =
yj V yj <Áy
где G - источниковые слагаемые, то это уравнение можно разделить на конвективную
df = -u df dt дх
и некоивективиую
df = G dt = G
части. При этом g будет удовлетворять уравнению
дд дд dG ди dt дх дх д дх' которое также расщепляется на конвективную и неконвективную части:
I + = ° <7>
дд dG ди dt дх дх
Решение конвективной части (7) осуществляется простым дифференцированием
f
* я f* dRliD[ Л^ Ci + C2Í + Сзе2 g =df =—(х - uAt) = 1 + ^—
Co + Cii + C2C2 + Сз£3
- aiPi-
(1 + aiM)2
Способ решения (8) обсуждается ниже при описании неконвективной стадии расчета.
В двумерном случае интерполяционная функция имеет вид [65]
Ч? (х,У) = I Е «р*вР,чXрУ" I £ Х 1ху 1у,
\0<р+д<1 у 0<1-х+1у<3
Ще X = X - ХН и У = у - уг,
« 0,0 = 1, в 0,0 = 1,
«1,0 = 1 при дхг,з дХгпР,о < 0, «1,0 = 0 при дхгпр] > 0,
«0,1 = 1 при дуг ^ дуг^пр < 0, «0,1 = 0 при дуг] дуЦпр > 0, в1,0 = [ |(Зх - д-хг,о )/(дхгпр,о - ^х)| - 1 ] А- 1,
в0,1 = [|(Зу - дуг] )/(дуг]пр - Зу)| - 1]А-1,
&0,0 = Д], С1,0 = дхг] + «1,0@1,01г,з, С0,1 = дуг] + «0,1в0,1/гл, С2,0 = [ (1 + «1,0в1,0Ах)/гпр] - С*0,0 - С^Ах ] А- 2 - Сз,0Ах, С0,2 = [ (1 + «0дв0дАу)1г]пр - С0,0 - С0дАу ] А- 2 - С0,зАу,
Сз,0 = [ (1 + а\,ов1,0^х)(дхгир,о — $х) + 9хг,о — ^^ Д- 2,
Со,3 = [ (1 + а0,1@0,1Ду )(дуг,оир — Яу) + дуг,2 — ] Д- 2,
, 1/гир, 3 + (1 + «1,0^1,0Дх)ду гир,3
+ [а1,0в1,0/г,3ир + (1 + ОЮдАцДу )дхгзир] Д-1 +
,0] г ,3ир I I '-•0,1^0,1" у л)хг,]ир ] "у
+ С0,0 — I /, ар,д вр,д /гир,]ир
\0<р+д<1 /
Д-1 Д-
+ Сз,0ДХД-1 + С0,зДу Д-1 + С2,0ДхД-1 + С0,2Ду Д-1,
С2,1 = [ 0-0,1^0,1/гир,з + (1 + а1,0в1,0Дх)дугир,] — С0>,1 ] Д-2 — С*1дД- 1,
С1,2 = [ &1,0@1,0 /*3ир + (1 + ^0,1^0,1 Ду)д-хг,зир — ^1,0] Д- 2 — С*1дД- 1.
Здесь гир = г — ^и > 0, тр = г + < 0, зир = 3 — ^и > 0,
Зир 3 + 1 При < Дх хгир,3 хг,3, Ду Уг,3ир У*,3, дхг,3 Вх/г,3,
дуг,з = ду/г,з, Ях = (/гир,з — /г,3)Д-\ Яу = (/г,зир — /г,3)Д^- На конвективной стадии значения искомой функции / и ее пространственных производных дх, ду определяются выражениями
/*,3 = (хг,з — иДЬ, уг,з — vДt),
вяз
д*хг,з = —вх~ (хг,з — uДt' Уг,з —
вя2П
Ву
д*уг,з = вгу (хг,3 — иД1,Уг,3 —
где и = (иг,3 + и*3)/2, V = + )/2, а и**, V** - скорости в точке (х*,3 — АЬ, Уг,3 — Д£). Значения и*г*3, предварительно рассчитываются с помощью Т1СЕР при и = и*3, V = . Для вычисления производных дх, ду на неконвективной стадии используется явная конечно-разностная аппроксимация, которая будет приведена ниже.
Для расчета функции-идентификатора у метод Т1С1Р применяется в сочетании со следующим приемом, позволяющим максимально сократить численную диффузию и немонотонность решения в окрестности больших градиентов. Вместо уравнения (2) решается уравнение
ВГ
ж + и ^ = °,
Г ^[(1 — е)п(у — 0.5)], (9)
параметр е определяет ширину размывания переходной области с ненулевым градиентом у. Функция-идентификатор у определяется через Г посредством обратного преобразования. Такой подход позволяет уменьшить переходную область между контактирующими средами до размеров одной-двух ячеек. Он является весьма эффективным в случае изменения связности контактных границ, в частности когда границы делятся, пересекаются или сливаются. Нужно отметить, что способ снижения численной диффузии (9) применим только в случае четко определенного диапазона изменения искомого параметра. Так, он может быть использо-
у
областью ее изменения [0,1], но не может применяться для скорости и, давления p, а в случае сжимаемой жидкости - и для плотности g.
На конвективной стадии метод RCIP применяется для решения уравнений переноса всех искомых параметров g, u, p, ^юш F (при использовании тангенциального преобразования (9)) и определения промежуточных значений пространственных производных от них. Окончательные значения производных рассчитываются на неконвективной стадии при решении набора уравнений типа (8). Конкретный вид каждого из этих уравнений определяется параметром, которому они соответствуют.
Для расчета неконвективной (акустической) части (4) применяется метод UP (Unified Procedure) [43]. Он является обобщением на сжимаемые среды методов расчета несжимаемых сред, использующих уравнение для давления (или для коррекции давления) таких, как MAC [66]. SMAC [67]. SIMPLE [68]. Одна из первых попыток такого обобщения предпринята в методе ICE [69]. Однако этот метод, как и его модификация PISO [70]. основан на консервативном подходе, при котором возникают упомянутые во введении трудности при расчете течений с контактными границами. В методе CIP-CUP (ClP-Combiried UP) применяются формулировка уравнений движения в неконсервативных переменных (плотность, скорость, давление) и расщепление расчета на конвективную и акустическую стадии. Это приводит к упрощению уравнения для давления, используемого на неконвективной стадии в методе UP. и открывает возможность расчета многофазных течений с большой разницей акустических имиедансов фаз.
Если давление представляется в виде p = p(g, s), то при малых изменениях плотности g и энтропии s для приращения давления Ap можно записать
-p = ( t ) Ag + ( dp 1 As.
в
Если же энтропия постоянна, то
Ap = ( ^ Ag = Of Ag. (10)
dg
На неконвективной стадии согласно (4)
Ag = gn+1 - g* = -Atg* V-un+1, (11)
Au = un+1 -u* = -At VP-. (12)
g* v ;
Применяя операцию дивергенции (V-) к уравнению (12) и выражая затем из него V • un+1 с учетом соотношения (10) и того, что Ap = pn+1 — p* , получаем из (11) следующее уравнение:
( Vpn+1 \_ pn+1 - p* , V • u* V Л g* ) = g*CS2At2 + At . (láj
Описанная процедура может применяться для расчета течений с широким диапазоном проявления сжимаемости среды. Кроме того, она обеспечивает непрерывность отношения Vp/g (которое в соответствии с (12) можно рассматривать как ускорение) на границе раздела сред даже при разнице плотностей порядка 103, в отличие от методов ISE и PISO, которые этого не гарантируют.
В настоящей работе используется модификация метода UP. предложенная в [71] с целыо повышения устойчивости и сходимости решения. Модификация основана
на схеме предиктор-корректор. На шаге предиктора находится промежуточное поле скорости по известному полю давления
й-и* = -Д* (14)
в*
Вычитая (14) из (12). получаем
и"+1 -Й = -д* ^ (15)
в*
где Sp = рп+1 -р* - коррекция давления. Применяя операцию дивергенции к уравнению (15) и выражая V - ип+1 из третьего уравнения из (4), приходим к следующему уравнению для коррекции давления:
fVSp\ 5р V и
В настоящей работе уравнение (16) решается методом последовательной верхней релаксации.
На шаге корректора применяются следующие соотношения:
вп+1 = в* - в* V - ип+1, ип+1 = й - д* + д*д„,
в*
рп+1 = р* + Sp + ДtQp.
В этих соотношениях слагаемые Д* Qu и Д* Qp отвечают за обеспечение искусственной вязкости, которая нужна при расчетах течений с сильными ударными волнами, поскольку схема С1Р-С11Р не является консервативной. Эти слагаемые определяются следующими выражениями [58]:
Qu = -1 Vqv, в
Qp = -(7 - l)qvV - и, qv = ^ в (-СзДи + 1тт1Ди2
Ди = шт(0, А V - и), А = Дх = Ду.
Коэффициент искусственной вязкости ^ обычно принимается для газа в интервале (0.5, 1), а для жидкости - в интервале (0, 0.5).
При обновлении пространственных производных от искомых параметров на неконвективной стадии согласно (8) сложностей с определением пространственных производных от источпиковых слагаемых О можно избежать. Действительно, пусть в одномерном случае на неконвективной стадии уравнение д//д* = О аппроксимируется выражением
/„+1 = /* + О.д*.
Тогда для дО/дх в уравнении (8) вместо аппроксимации
°г+1 - Ог-1
2Дх,-
можно использовать
хп+1 Г* Л*п+1 , Г*
^¿+1 ^¿+1 ^¿-1 + ^¿-1
2ДíДx¿
С учетом этого аппроксимацию уравнения (8) можно записать как
п+1 * -/*П+1 Г* ХП+1 , Г* / \
= /¿+1 - /¿+1 - /¿-1 + /¿-1 - * К+1 - ^¿-1) Д; 2ДtДx¿ 0 2Дx¿ '
В двумерном случае для обновления пространственных производных =
= / на неконвективной стадии используется следующая конечно-разностная аппроксимация [72]:
/п+1 _ /* _ /п+1 + /*
п+1 = * + ^ ¿+1,3 ^ ¿+1,3 ^¿-1,3 + ^¿-1,3
= ^¿,з + 2Д^,3
* (^¿+1,3 _ ^¿-1,3 )Д* _ * (^¿+1,3 _ ^¿-1,3 )Д* 2Дx¿,д• ,
/п+1 _ /* _ /п+1 + /*
п+1 _ Л* | ^¿,3+1 ^¿,3+1 ^¿,3-1 + ^¿,3-1
(17)
Оп+1 = о* . . +
^¿,з ^¿,з + 2Дy¿,j
_ * (^¿,3+1 _ ^¿,3-1)Д£ _ * (^¿,3+1 _ ^¿,3-1)Д^
^ 2Дy¿,j 2Дуг,3 '
В дальнейшем сочетание метода С1Р-С11Р с вышеописанной методикой выявления контактной границы на эйлеровых сетках будет упоминаться как С1Р-С11Р1.
2. Верификация методики расчета
Проверка правильности работы алгоритма и программы расчета осуществлялась на примере задачи о вращении твердого диска с вырезом [73] и о распаде разрыва на контактной границе между газом и жидкостью.
2.1. Вращение твердого диска с вырезом. Конвективная часть уравнений движения сред (3) вместе с уравнением для функции-идентификатора (2) представляют собой набор уравнений переноса для параметров и, р, у. Тестирование реализации метода ИСТР. применяемого для их численного решения, осуществляется с использованием задачи о вращении твердого диска с вырезом в области 0 < х < 100, 0 < у < 100 вокруг точки х = 50, у = 50 [73].
Диск вращается вместе со средой, безразмерное поле скорости которой имеет вид и = (у _ 50), V = _ (х _ 50). Рассматривается один оборот вокруг центра
вращения. Изменение положения диска определяется решением уравнения пере-
у
у
перепоса для у решается та равномерной сетке, со стоящей из 100 х 100 ячеек с пространственными шагами Дх = Ду = 1. Шаг то времени Д; = 0.01. Результаты расчетов приведены на рис. 2.
Положение диска определяется по пространственному полю функции-идонти-фикатора у: диск находится в области, где у = 1, а граница диска - там, где
0 < у < 1. Согласно начальной дискретизации у изменяется от 0 до 1 в области
у
диска выглядит как четкий контур (рис. 2, о). При решении уравнения переноса (2) методом Т1С1Р с малой, но все же ненулевой численной диффузией изначально
Рис. 2. Изолинии функции описывающие поверхность твердого диска с вырезом, в начале его вращения (а) и после завершения одного оборота, рассчитанные методами RCIP (б), RCIP с тангенциальным преобразованием (9) (е), с использованием схем против потока (г) и Лакса-Вендроффа (à) с тангенциальным преобразованием. Результаты (г) и (д) получены в [74]. В начальный момент времени (а), для которого представлена также часть расчетной области с полем скорости, все изолинии слились в одну жирную кривую, представляющую собой начальную границу диска
резкий профиль ^ на границе диска размывается, так что неопределенность положения границы в конце одного оборота увеличивается (рис. 2, б). Однако монотонность решения сохраняется, и форма границы, несмотря на размывание, имеет вид, подобный начальному. Если метод RCIP применяется в сочетании с антидиффузионным тангенциальным преобразованием (9), то профиль диска в поле изолиний функции ^ в конце вращения (рис. 2, в) практически полностью совпадает с начальным. Поэтому можно заключить, что используемый в настоящей работе метод расчета конвективных слагаемых позволяет получать монотонные решения в области больших градиентов и характеризуется малой численной диффузией. При расчете конвективной части уравнений движения сред метод RCIP применяется для вычисления параметров u, p. При дополнительном использовании тангенциального преобразования для вычисления функциии-идентификатора ^ граница между средами на протяжении всего расчета размывается не более чем на одну-две ячейки. Нужно отметить, что применение тангенциального преобразования улучшает не все схемы [74]. Так, и в сочетании с этим преобразованием схема "против потока" (рис. 2, г) и схема Лакса-Вендроффа (рис. 2, д) дают очень большие искажения границы диска. При этом схема "против потока" остается диффузионной, а схема Лакса-Вендроффа - немонотонной.
2.2. Распад разрыва на границе газ— жидкость. Для тестирования полного алгоритма расчета, состоящего из конвективной и акустической стадии, используется одномерная задача о распаде плоского разрыва со скачком давления на границе газ-жидкость.
Расчет проводится в области [0, 5х0] ? где х0 - начальное положение разрыва, на равномерной сетке из 600 ячеек с шагом Ах = 0.03 м и числом Куранта 0.8. Параметры сред p, u, y, B) по разные стороны разрыва заданы следующим образом (1 кг/м3; 1 атм; 0; 1.4; 0) для газа при х < х0 и (1.2• 103кг/м3; 5• 103 атм; 0; 7; 3 • 103 атм) для жидкости при х > хо . Коэффициенты искусственной вязкости полагаются равными 0.6 для газа и 0 - для жидкости. Результаты расчетов и аналитическое решение представлены на рис. 3.
Рис. 3. Распад разрыва на границе газ (слева) - жидкость (справа); хо - начальное положение разрыва (отмечено вертикальной штриховой лилией). Сплошные кривые аналитическое решение Римана, символы • - численное решение методом С1Р-СиР1, штриховые кривые численное решение методом С1Р-СиР1 без искусственной вязкости. Цифрами отмечены: 1 - движущаяся влево контактная граница; 2 - волна разрежения, распространяющаяся по жидкости вправо; 3 - ударная волна, распространяющаяся по газу влево
Как видно, при численном решении задачи методом С1Р-С11Р1 без использования искусственной вязкости положение контактной границы (кривые 1) и волны разрежения в жидкости (кривые 2) хорошо согласуется с аналитическим решением Римана, а в окрестности фронта ударной волны в газе (кривые 3) возникают осцилляции давления и положение ее фронта получается неверным (штриховые кривые). Корректный расчет ударной волны удается реализовать только с использованием искусственной вязкости (символы •).
3. Иллюстрация возможностей метода
Для иллюстрации возможностей реализованной методики расчета рассматривается задача об ударном воздействии высокоскоростной струи жидкости по твердой стенке. Предполагается, что параметры струи соответствуют параметрам струйки жидкости, возникающей на поверхности кавитационного пузырька при его схлопы-вании либо в контакте со стенкой, либо на небольшом удалении от нее. Эта задача связана с проблемой кавитационного разрушения поверхностей тел, которая представляет значительный интерес в связи со многими практическими приложениями. Несмотря на большое количество исследований воздействия кавитационных пузырьков на стенку [40, 44, 75 77], ряд важных аспектов такого воздействия до сих пор остается неизученным. К ним, в частности, относится рассматриваемый в настоящей работе режим ударного воздействия, при котором в финальной стадии сжатия пузырька в жидкости возникают ударные волны.
Как показывают эксперименты [78 81], в процессе несферичного охлопывания пузырька у стопки возможен сценарий о формированием направленной к стенке высокоскоростной струйки жидкости, которая пересекает полость пузырька и бьет
я
>>
а н
Я
V
V
////////////////////////////////////)//)///&//////?/?//////////////
Рис. 4. Схема ударного воздействия пузырька па стопку
либо по противоположной частн поверхности пузырька (еслн пузырек не примыкает к стенке), либо непосредственно по стенке (в случае примыкания). Расчеты [82] показывают, что наибольшие скорости струйки достигаются в случае, когда диаметр струйки мал по сравнению с размерами пузырька (рис. 4. а). В связи с этим на начальном этапе исследования ударного воздействия пузырька на твердую поверхность влиянием верхней и боковых границ пузырька можно пренебречь. В результате ударное воздействие струйки можно моделировать как удар столбика жидкости по стенке (рис. 4. б, в) или по жидкой прослойке на ней (рис. 4. г).
В настоящей работе рассматривается осесимметричиая задача об ортогональном ударе цилиндрической струи жидкости с плоским (рис. 4. б) или полусферическим (рис. 4. е) концом по неподвижной плоской твердой стенке или по тонкой жидкой прослойке на ней (рис. 4. г). Струя считается бесконечно длинной и окруженной неограниченным объемом газа. Эффекты вязкости, теплопроводности жидкости и газа, а также влияние поверхностного натяжения не учитываются. Для жидкости начальные значения параметров полагаются следующими: д = 103 кг/м3, 7 = 7.15, В = 3047 атм, для газа: д = 1 кг/м3, 7 = 1.4, В = 0. В начальный момент времени давление всюду равно 1 атм. Радиус струи 25 мкм, 500 /
кумулятивной струйки жидкости, образующейся на поверхности кавитационного пузырька в воде при комнатных условиях в том случае, когда пузырек примыкает к телу и в начале охлопывания имеет вид сплюснутого сфероида (эллипсоида вращения с отношением полуосей ~ 0.84). Коэффициент искусственной вязкости еу в жидкости полагается равным 0.4, а в газе — 0.9. Применяются равномерные или сгущающиеся у стенки расчетные сетки с соотношением числа ячеек вдоль координатных направлений примерно равным 400 х 300. Число Куранта варьировалось в интервале 0.8 ^ 0.9.
3.1. Удар струи жидкости с плоским концом по твердой стенке.
Схема рассматриваемого удара струи жидкости с плоским концом по твердой стенке изображена на рис. 4, б. На рис. 5 представлено сопоставление расчетных сеток настоящей работы (а) и работы [77] (б), а также формы струи в некоторый момент времени на стадии ее растекания. В работе [77] сетка является адаптивно-подвижной, а влияние окружающего газа не учитывается. Для наглядности сетки разрежены в 8 раз по каждому направлению.
В настоящей работе граница струи определяется неявно, ее положение, рассчитанное на эйлеровой сетке, показано жирной кривой на рис. 5, о. Эта граница образована набором близкорасположенных изолиний функции-идентификатора у, подобно тому, как кривая ударной волны образуется близкорасположенными изолиниями давления. В работе [77] расчетная сетка подстраивается под границу струи, которая всегда является координатной линией сетки (рис. 5, б), а уравнения движения среды решаются с применением модификации метода Годунова второго порядка точности.
0.6-
г/Я
а
0.6-в=
0.6 г/Я 1.2
0.6 г/Я 1.2
Гис. 5. Гасчетпые сетки в области взаимодействия струи с плоским концом со стенкой и граница струн (жирная кривая) па стадии растекания: а эйлерова сетка и рассчитанная па пей методом С1Р-СиР1 граница струи: б адаптивно-подвижная сетка работы [77], правой границей которой является граница струи. Левая граница сеток является осыо симметрии
Рис. 6. Прореженное поле скорости, изолинии давления и изолинии фупкции-идептифи-катора (а) в момент времени, которому соответствует рис. 5: а расчет методом С1Р-СиР1: б результаты из работы [77]. Цифрами отмечены: 1 ударная волна: 2 волна разрежения: 3 граница струи
Рис. 7. Распределения давления (а) и горизонтальной составляющей скорости (б) па стенке в последовательные моменты времени ¿1, , ¿з в случае удара струи с плоским концом: символы • - расчет методом С1Р-С11Р1; сплошные кривые — результаты из работы [77]
На рис. 6 представлено сопоставление изолиний давления и полей скорости, полученных в настоящей работе методом С1Р-С11Р1 (о) и в работе [77] (б) в тот же момент времени, что и на рис. 5. На рис. 7 показано сравнение полученных методом С1Р-С11Р1 н в [77] распределений давления и скорости на стенке для нескольких моментов времени. Следует отметить, что ударная волна 1 (рис. 6, о), возникающая в струе в момент удара по стенке, является довольно интенсивной. Скорость ее распространения составляет ~ 1900 м/с (для сравнения скорость звука в невозмущенной струе ~ 1500 м/с), а перепад давления па ее фронте ~ 104 атм. Согласно рис. 6, 7 применение искусственной вязкости позволяет корректно определить как положение ударной волны, так и значение давления за ее фронтом. Таким образом,
Рис. 8. Расчетные сетки в области взаимодействия струн с полусферическим концом со стенкой в конце стадии патекапия струи: а сетка настоящей работы: б адаптивпо-по-движиая сетка работы [77]. Жирная кривая граница струи
Рис. 9. Прореженное поле скорости, изолинии давления и изолинии фупкции-идептифи-катора (а) в конце стадии патекапия струи в момент времени, которому соответствует рис. 8: а расчет методом С1Р-СиР1: б результаты из работы [77]. Цифрами отмечены: 1 ударная волна: 2 волна разрежения: 3 граница струи
несмотря на существенное различие методик работы [77] и настоящей работы, согласование результатов вполне удовлетворительное, что свидетельствует о работоспособности настоящей реализации метода С1Р-С11Р1.
3.2. Удар струи жидкости с полусферическим концом по твердой стенке. Следующим шагом в моделировании ударного воздействия на стенку кумулятивной струйки, образующейся при несферическом схлопывании пузырька, является решение задачи об ударе цилиндрической струи с полусферическим концом по твердой стенке. Эта задача ближе к реальным условиям по сравнению с рассмотренным выше ударом цилиндрической струи с плоским концом, но и более сложна для численного моделирования (рис. 4. в). Одна из особенностей этой задачи заключается в том. что начальное натекание струи сменяется в некоторый момент ее растеканием с образованием тонкой пристеночной струйки [12. 77]. При этом использование адаптивно-подвижных сеток [77] усложняется в связи с их сильными деформациями, необходимостью редискретизации расчетной области с применением сеток различающейся топологии и интерполяции решения с одной сетки на другую.
Результаты расчетов удара струи жидкости с полусферическим концом по твердой стенке приведены на рис. 8 11. На рис. 8 показано сопоставление расчетных сеток, применяемых в настоящей работе (о) ив расчетах [77] (б), а также формы струи в конце стадии ее натекания. за которой следует растекание. На рис. 9 дается соответствующее рис. 8 сопоставление изолиний давления и полей скорости в конце стадии натекания струи на стенку. Как видно, положение границы струи (кривая 3), а также положение и криволинейная форма ударной волны в струе (кривая 1), полученные с помощью метода С1Р-С11Р1, хорошо согласуются с результатами [77]. При этом важно отметить, что в настоящей работе начальный
Гис. 10. Гастекание струи с полусферическим концом: 1 - ударная волна; 2 - волна разрежения; 3 - граница струи
Гис. 11. Радиальные распределения давления (в) и горизонтальной составляющей скорости (б) на стенке в последовательные моменты времени Ь1,Ь2,Ь3,Ь4 в случае удара струи с полусферическим концом: символы • - расчет методом С1Р-СиР1; сплошные кривые — результаты из работы [77]. Моменту £3 в конце стадии натекания струи соответствуют рис 8, 9, моменту Ь4 на стадии растекания соответствует рис. 10
полукруглый профиль струи аппроксимируется на эйлеровой сетке отрезками прямых. что вносит определенные погрешности в результаты. Для уменьшения этих погрешностей необходимо измельчать сетку. В настоящей работе с целыо сокращения объема вычислений осуществлялось сгущение расчетной сетки у стенки.
В работе [77] была рассчитана только стадия натекания струи. Расчет ее дальнейшего растекания не осуществлялся в связи со значительными сложностями перестройки сеток. В то же время метод С1Р-С11Р1 позволяет рассчитывать не только всю стадию натекания струи, но и ее последующее растекание (с формированием тонкой пристеночной струйки) (рис. 10). При этом каких-либо дополнительных трудностей не возникает. На рис. 11. а приведено сравнение полученных в настоящей работе методом С1Р-С11Р1 и в работе [77] пристеночных распределений давления. Представлены четыре последовательных момента времени, включая момент растекания струи (которому соответствует рис. 10). рассчитанный только методом С1Р-С11Р1. Можно отметить хорошее согласование с результатами из работы [77] как значений давления за фронтом ударной волны (кривые, соответствующие моменту ¿1), так и распределений давления на стенке в моменты ¿1, ¿2, ¿3- Некоторая немонотонность кривых давления, полученных с использованием С1Р-С11Р1. а также меньшее по сравнению с [77] максимальное значение давления на пике в окрестности контакта границы струи со стенкой в момент 13 могут быть обусловлены несовпадением границы струи с узлами сетки и аппроксимацией этой границы отрезками прямых линий. При измельчении сетки эти отличия от результатов работы [77] уменьшаются. На рис. 11. б сравниваются профили горизонтальной составляющей скорости на стенке в конце стадии натекания струи (для большей
1
1
а
б
0.5-
0
1
3
о
о
0.5
1 т 1-5
о
0.5
1 г/Я 1.5
Рис. 12. Удар струи с плоским концом по топкому слою жидкости па стопке. Изолинии давления и фупкции-идептификатора: а па начальной стадии удара: б па стадии развития выплеска с отрывом маленькой капли. Цифрами отмечены: 1 ударная волна: 2 волна разрежения: 3 границы струи и жидкой прослойки с выплеском: 4 капля
наглядности только в момент, которому соответствуют рис. 8. 9). полученные в настоящей работе методом С1Р-С11Р1 и в работе [77]. Наряду с хорошим согласованием положения максимума горизонтальной скорости с результатами [77]. можно отметить немонотонный характер кривой настоящей работы и большее значение максимума. Эти отличия также могут быть обусловлены особенностями аппроксимации криволинейной границы контакта на прямолинейной сетке и уменьшаются по мере ее измельчения.
Из рис. 8. 9 и 11 можно заключить, что. несмотря на погрешности, связанные с аппроксимацией гладких кривых па эйлеровой сетке отрезками прямых линий, согласование результатов расчетов методом С1Р-С11Р1 с расчетами работы [77] вполне удовлетворительное. При этом следует подчеркнуть, что метод С1Р-С11Р1 позволяет избежать характерных для метода адаптивно-подвижных сеток трудностей. связанных с большим деформированием сеток.
3.3. Удар струи жидкости с плоским концом по жидкой прослойке на твердой стенке. Когда пузырек не примыкает к стенке, а находится на относительно небольшом (порядка его радиуса) удалении от нее. струйка, формирующаяся при несферическом схлопывании пузырька, ударяет по жидкой прослойке между пузырьком и стенкой. При этом воздействие на стенку осуществляется посредством ударных волн, возникающих в жидкой прослойке в момент удара. Этот вариант для струи с плоским концом схематично показан на рис. 4. г. Отличие от постановки рассмотренной в п. 3.1. задачи об ударе струи с плоским концом по стенке состоит лишь в том. что в данном случае на стенке имеется тонкий слой жидкости, которая идентична жидкости в струе. В то же время из-за сильного деформирования контактной границы данный случай гораздо сложнее для численного моделирования на основе адаптивно-подвижных сеток. Поэтому эта задача не была решена авторами работы [77].
На рис. 12 представлены результаты расчета с использованием метода С1Р-С11Р1 удара струи жидкости с плоским концом по жидкой прослойке на твердой стенке. Как видно, данный случай характеризуется не только сильной деформацией границы контакта газ жидкость при формировании выплеска жидкости в результате удара струи по жидкой прослойке (рис. 12. о), но и изменением связности контактной границы при срыве маленькой капли жидкости с вершины выплеска (рис. 12. б). Метод С1Р-С11Р1 позволяет рассчитывать такие межфазные взаимодействия естественным образом, не выходя за рамки имеющихся численных алгоритмов.
Заключение
Представлен краткий обзор работ по численному моделированию задач контактного взаимодействия сжимаемых сред, в том числе с большой разницей акустических имподансов при сильных деформациях контактной границы и при наличии в средах ударных волн. Приведены основные положения реализованной методики расчета таких задач на основе метода CIP-CUP, в которой применяются эйлеровы сотки без явного выделения контактных границ. Такой подход позволяет избежать проблем с построением расчетных сеток как при сильных деформациях границ, так и при изменении их топологии. Представлены результаты расчетов модельных задач, имеющих аналитическое решение, подтверждающие эффективность созданного алгоритма. Возможности реализованной методики проиллюстрированы на задачах об ударе осеспмметрпчной высокоскоростной струи жидкости по жесткой стенке и по топкому слою жидкости на жесткой стенке. Полученные результаты сопоставлены с известными численными решениями, рассчитанными методом адаптивно-подвижных сеток с явным выделением межфазной границы. В некоторых случаях сопоставление проведено лишь на начальной стадии удара, поскольку дальнейшее применение метода адаптивно-подвижных сеток становилось проблематичным. Показано удовлетворительное согласование результатов.
Работа выполнена в рамках программы РАН. гранта президента РФ (МК-2712.2011.1) и при поддержке РФФИ (проект Л* 12-01-00341-а).
Summary
A.A. Aganin, T.S. Guseva. Numerical Simulation of Contact. Interaction of Compressible Fluids 011 Eulerian Grids.
A short, review 011 numerical simulation of the problems of contact, interaction between compressible media with a strong difference in their acoustic impedances, with large deformations of contact, boundary and with shock waves present, in t.hem is given. The main ideas underlying a numerical algorithm developed for solving such problems 011 the basis of the known CIP-CUP (Constrained Interpolation Profile Combined Unified Procedure) approach using Eulerian grids without, explicit, separation of contact, interfaces, and governing equations of compressible fluid flow in terms of nonconservat.ive variables are described. Robustness of the developed algorithm is demonstrated by the results of computation of some test, problems having analytical solutions. The capabilities of the realized technique are illustrated by its application to the problems of impact, of an axially-symmet.ric high-speed liquid jet. 011 a rigid wall and 011 a thin liquid layer 011 a rigid wall. The obtained results are in satisfactory agreement, with the known numerical solutions computed by the method of adapt.ively moving meshes with an explicit, separation of the interphase boundary.
Key words: multiphase compressible flows, large deformation of interface, Eulerian grid, interface capturing, CIP-CUP method.
Литература
1. Benson D.J. Computational methods in Lagrangian and Eulerian liydrocodes // Comput.. Method. Appl. M. 1992. V. 99, No 2 3. P. 235 394.
2. Mair H. U. Review: liydrocodes for structural response to underwater explosions // Shock Vib. 1999. V. 6, No 2. P. 81 96.
3. Allaire G., Clere S., Kukh S. A five-equation model for the simulation of interfaces between compressible fluids // J. Comput.. Pliys. 2002. V. 181, No 2. P. 577 616.
4. Бураго Н.Г., Кукудж.аиов B.H. Обзор контактных алгоритмов // Изв. РАН. МТТ. 2005. 1. С. 45 87.
5. Kadioglu S.Y., Sussman М., Osker S., Wright J.P., Капу M. A second order primitive preconditioner for solving all speed multi-pliase Hows // J. Comput. Pliys. 2005. V. 209, No 2. P. 477 503.
6. Garimella R., Dyadeehko V., Swartz В., Shashkov M. Interface reconstruction in multi-fluid, multi-pliase How simulations // Proc. 14t.li Int.. Meshing Roundt.able. Springer. 2005. P. 19 32.
7. Johnsen E., Culunius Т. Implementation of WENO schemes in compressible mult.icomponent. flow problems // J. Comput. Pliys. 2006. V. 219, No 2. P. 715 732.
8. Hu X.Y., Khuu B.C., Adams N.A., Huang F.L. A conservative interface method for compressible flows // J. Comput. Pliys. 2006. V. 219, No 2. P. 553 578.
9. Kim J., Lowengrub J. Interfaces and mult.icompoiieiit. fluids // Encyclopedia of Mathematical Physics. Elsevier, 2006. P. 135 144.
10. Dunghung W., Ning Z., Ou H., Jianming L. A Ghost Fluid based Front. Tracking Method for mult.imedium compressible flows // Acta Math. Sei. 2009. V. 29, No 6. P. 1629 1646.
11. So K.K., Hu X.Y., Adams N.A. Anti-diffusion method for interface steepening in two-phase incompressible flow // J. Comput.. Pliys. 2011. V. 230, No 13. P. 5155 5177.
12. Чнжов A.B., Шлтд'т A.A. Высокоскоростной удар капли о преграду // Журп. теорет. физики. 2000. Т. 70, Вып. 12. С. 18 27.
13. Головачев Ю.П., Ноткииа Е.А., Чижов A.B., Шлтд'т A.A. Расчет ударпо-волповых течений со свободными поверхностями // Журп. вычисл. матем. и матем. физики. 2001. Т. 41, Л» 1. С. 157 167.
14. Годунов С.К., Забродин A.B. и др. Численное решение многомерных задач газовой динамики. М.: Наука, 1976. 400 с.
15. Unverdi S.O., Tryggvasun G. A front-tracking method for viscous, incompressible, multi-fluid flows // J. Comput.. Pliys. 1992. V. 100, No 1. P. 25 37.
16. Cueehi J.-P., Saurel R. A Riemann problem based method for the resolution of compressible mult.imat.erial flows // J. Comput.. Pliys. 1997. V. 137, No 2. P. 265 298.
17. Glimm J., Grove J.W., Li X.-L., Shyue K.-M., Zhang Q., Zeng Y. Three dimensional Front. Tracking // SIAM J. Sei. Comput.. 1998. V. 19, No 3. P. 703 727.
18. Glimm J., Grove J., Li X.-L., Tan D. Robust, computational algorithms for dynamic interface tracking in three dimensions // SIAM J. Sei. Comput.. 2000. V. 21, No 6. P. 2240 2256.
19. Haller К.К., Ventikos Y., Poulikakos D. Computational study of high-speed liquid droplet, impact. // J. Appl. Pliys. 2002. V. 95, No 2. P. 2821 2828.
20. Terashima H., Tryggvason G. A front-t.racking/ghost.-fluid method for fluid interfaces in compressible flows // J. Comput.. Pliys. 2009. V. 228, No 11. P. 4012 4037.
21. Harlow F., Amsden A. Fluid dynamics: a LASL monograph. Springfield, Virginia: Nat.. Teclm. Inform. Serv., 1971. 115 p.
22. Abgrall R., Kami S. Computations of compressible mult.ifiuids // J. Comput.. Pliys. 2001. V. 169, No 2. P. 594 623.
23. Hirt C.W., Niehols B.D. Volume of fluid (VOF) method for the dynamics of free boundaries // J. Comput.. Pliys. 1981. V. 39, No 1. P. 201 225.
24. Wemmenhove R., Veldman A.E.P., Luppes R., Bunnik T. Application of a VOF method to model compressible t.wo-pliase flow in sloshing tanks // Proc. of OMAE2008. 2008. P. 603 612. Paper Л» OMAE2008-57254.
25. Петров Н.В. Исследование локального эперговыделепия в воде вблизи свободной поверхности // Науч.-техп. ведом. СПбГПУ. 2012. Л' 1. С. 92 96.
26. Rider W.J., Kothe D. Reconstructing volume tracking // J. Comput. Pliys. 1998. V. 141, No 2. P. 112 152.
27. Renardy Y., Renardy M. PROST: a parabolic reconstruction of surface tension for the volume-of-fluid method // J. Comput. Pliys. 2002. V. 183, No 2. P. 400 421.
28. Ubbink O., Issa R.I. A method for capturing sharp fluid interfaces on arbitrary meshes // J. Comput. Pliys. 1999. V. 153, No 1. P. 26 50.
29. Cassidy D.A., Edwards J.R., Tian M. An investigation of interface-sharpening schemes for multi-phase mixture flows//J. Comput. Phys. 2009. V. 228, No 16. P. 5628 5649.
30. Osher S., Sethian J. Front propagating with curvature dependent speed: Algorithms based on Hamilton Jacobi formulations // J. Comput. Pliys. 1988. V. 78, No 2. P. 12 49.
31. Sethian J.A., Smereka P. Level Set Methods for fluid interfaces // Annu. Rev. Fluid Mech. 2003. V. 35. P. 341 372.
32. Kami S. Multicomponent. flow calculation by a consistent primitive algorithm // J. Comput. Pliys. 1994. V. 112, No 1. P. 31 43.
33. Гапоненко Ю.А. Численное моделирование газовой кумуляции продуктов взрыва при детонации плоских параллельных зарядов // Вычисл. технологии. 2000. Т. 5, Л» 4. С. 31 39.
34. Abgrall R. How to prevent pressure oscillations in multicomponent flow calculations: A quasi-conservative approach // J. Comput. Pliys. 1996. V. 125, No 1. P. 150 160.
35. Anderson D.M., MeFadden G.B., Wheeler A.A. Diffuse-interface methods in fluid mechanics // Aim. Rev. Fluid Mech. 1998. V. 30. P. 139 165.
36. Lowengrub J., Truskinuvsky L. Quasi-incompressible Calm-Hilliard fluids and topological transitions // Proc. Roy. Soc. Loud. A. 1998. V. 454. P. 2617 2654.
37. Larrouturou B. How to preserve the mass fraction posit.ivit.y when computing compressible multi-component, flows // J. Comput. Pliys. 1990. V. 95, No 1. P. 59 84.
38. Shyue K.M. An efficient, shock-capturing algorithm for compressible multicomponent. problems // J. Comput.. Pliys. 1998. V. 142, No 1. P. 208 242.
39. Abgrall R., Saurel R. Discrete equations for physical and numerical compressible multiphase mixtures // J. Comput.. Pliys. 2003. V. 186, No 2. P. 361 396.
40. О sternum A., Dular M., Siruk B. Numerical simulation of a near-wall bubble collapse in an ultrasonic field // J. Fluid Sci. Teclm. 2009. V. 4, No 1 P. 210 221.
41. Shyue K.M. An anti-diffusion based Eulerian interface-sharpening algorithm for compressible t.wo-pliase flow with cavitation // Proc. 8t.li Int.. Symposium on Cavitation. Singapore, 2012. - Abstr. No 198. - URL: http://www.math.ntu.edu.tw/~shyue/ mypapers/kmsliyue_cav2012.pdf, свободный.
42. Ballil A., Julgam S., Nuwakuwski A.F., Nieuleau F.C.G.A. Numerical simulation of compressible t.wo-pliase flows using an Eulerian type reduced model // Proc. World Congress on Engineering, WCE 2012. London, 2012. P. 1835 1840.
43. Yabe Т., Xiao F., Vtsumi T. The constrained interpolation profile method for multiphase analysis // J. Comput.. Pliys. 2001. V. 169, No 2. P. 556 593.
44. Duihara R., Takahashi K. Numerical calculation of laser-produced bubble near a solid boundary until the second collapse // JSME. Int.. J. Ser. B. 2001. V. 44, No 2. P. 238 246.
45. Kawasaki К. Numerical model of 2-D multiphase flow with solid-liquid-gas interaction // Int. J. Offshore Polar Eng. 2005. V. 15, No 3. P. 198 203.
46. Kukh S., Allaire G. Numerical simulation of 2D t.wo-pliase flows with interface // Godunov Methods: Theory and Applications. 2001. P. 513 518.
47. Xiao F., Honma Y., Kono T. A simple algebraic interface capturing scheme using hyperbolic tangent function // Int. J. Numer. Methods Fluids. 2005. V. 48, No 9. P. 1023 1040.
48. Sun Y., Beckermann C. Sharp interface tracking using the phase-field equation // J. Comput. Phys. 2007. V. 220, No 2. P. 626 653.
49. Yabe Т., Wang P.Y. Unified numerical procedure for compressible and incompressible fluid // J. Phys. Soc. Japan. 1991. V. 60, No 7. P. 2105 2108.
50. Falkiw R., Liu X.-D., Osher S. A general technique for eliminating spurious oscillations in conservative schemes for multiphase and multispecies Euler equations // Int. J. Nonlinear Sci. Numer. Sim. 2002. V. 3, No 2. P. 99 106.
51. Kami S. Hybrid multifluid algorithms // SIAM J. Sci. Comput. 1996. V. 17, No 5. P. 1019 1039.
52. Glimm J., Li X.L., Liu Y.-J., Xu Z.-L., Zhao N. Conservative front tracking with improved accuracy // SIAM J. Num. Analysis. 2003. V. 41, No 5. P. 1926 1947.
53. Falkiw R., A slam Т., Merriman В., Osher S. A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the Ghost Gluid Method) // J. Comput. Phys. 1999. V. 152, No 2. P. 457 492.
54. Falkiw R. Coupling an Eulerian fluid calculation to a Lagrangian solid calculation with the ghost fluid method // J. Comput. Phys. 2002. V. 175, No 1. P. 200 224.
55. Nouryaliev R.R., Dinh T.N., Theofanous T.G. Direct numerical simulation of compressible multiphase flows: interaction of shock waves with dispersed multimaterial media // Proc. 5t.li Int. Conf. Multiphase Flow, ICMF!04. Yokohama, Japan, 2004. Paper No 494. URL: http://wnAfw.crss.ucsb.edu/music/LEVELO/ConferencesOpen/ Conferences.2004/DNS-ICMF04.pdf, свободный
56. Takahira H., Yuasa S. Numerical simulations of shock-bubble interactions using an improved Ghost Fluid Method // ASME 2005 Fluids Engineering Division Summer Meeting (FEDSM2005). 2005. No FEDSM 2005-77119. P. 777 785.
57. Nouryaliev R., Kadioglu S., Mousseau V. Fully-implicit, interface tracking for all-speed multifluid flows // Computational Fluid Dynamics 2008. Berlin: Heidelberg: Springer, 2009. P. 551 557.
58. Ogata Y., Yabe T. Shock capturing with improved numerical viscosity in primitive Euler representation // Comput. Phys. Commun. 1999. V. 119, No 2 3. P. 179 193.
59. Peng G., Ishizuka M., Hayama S. An improved CIP-CUP method for submerged water jet flow simulation // JSME Int. J. Ser. B. 2001. V. 44, No 4. P. 497 504.
60. Ida M. An improved numerical solver for the unified solution of compressible and incompressible fluids involving free surfaces. II. Multi-Time-Step integration and applications // Comput. Phys. Commun. 2003. V. 150, No 3. P. 300 322.
61. Nomura S., Nishida K. Numerical simulation of a single bubble rising in an ultrasonic standing wave field // Jpn. J. Appl. Phys. 2003. V. 42. P. 2975 2980.
62. Tong M., Browne D.J. Modelling compressible gas flow near the nozzle of a gas atomiser using a new unified model//Computers and Fluids. 2009. V. 38, No 6. P. 1183 1190.
63. Тапака N.. Maseguchi R., Ogawara Т. Improvement of conservative property and interlace resolution of mesli-based t.wo-pliase flow simulation algorithms for splashing fluid behavior // Nuc. Eng. Design. 2010. V. 240, No 12. P. 3984 3990.
64. Lee S.-J., Cho В.-G., Lee I. Two-dimensional unsteady aerodynamics analysis based on a multiphase perspective // Comput. Fluids. 2012. V. 53. P. 105 116.
65. Xiao F. A class of single-cell high-order semi-Lagrangian advect.ion schemes // Month. Weat.li. Rev. 2000. V. 128, No 4. P. 1165 1176.
66. Harlow F.H., Welch J.E. Numerical calculation of time-dependent, viscous incompressible flow with free surface // Pliys. Fluids. 1965. V. 8. P. 2182 2189.
67. Amsden A.A., Harlow F.H. A simplified MAC technique for incompressible fluid flow calculations // J. Comput.. Pliys. 1970. V. 6, No 2. P. 322 325.
68. Patankar S.V., Spalding D.B. A calculation procedure for heat., mass and momentum transfer in three-dimensional parabolic flows // J. Heat. Mass Transfer. 1972. V. 15. P. 1787 1806.
69. Harlow F.H., Amsden A,A, Numerical simulation of almost, incompressible flow // J. Comput.. Pliys. 1968. V. 3, No 1. P. 80 93.
70. Issa R.I. Solution of the implicitly discret.ised fluid flow equations by operator-splitting // J. Comput.. Pliys. 1986. V. 62, No 1. P. 40 65.
71. Yoon S.Y., Yabe T. The unified simulation for incompressible and compressible flow by the predictor-corrector scheme based on the CIP method // Comput.. Pliys. Commun. 1999. V. 119, No 2 3. P. 149 158.
72. Yabe T., Ishikawa T., Wang P.Y., Aoki T., Kadota Y., Ikeda F. A universal solver for hyperbolic equations by cubic-polynomial interpolation. II. Two- and three-dimensional solvers // Comput.. Pliys. Commun. 1991. V. 66, No 2 3. P. 233 242.
73. Zalesak S.T. Fully multidimensional flux-corrected transport, algorithms for fluids // J. Comput.. Pliys. 1979. V. 31, No 3. P. 335 362.
74. Yabe T., Xiao F. Description of complex and sharp interface with fixed grids in incompressible and compressible fluid//Comput. Math. Appl. 1995. V. 29, No 1. P. 15 25.
75. Blake J.R., Gibson B.C. Cavitation bubbles near boundaries // Annu. Rev. Fluid Mecli. 1987. V. 19. P. 99 123.
76. Аганин А.А., Халитоаа Т.Ф., Хисматуллииа, H.А. Метод численного решения задач сильного сжатия песферического кавитациошюго пузырька // Вычисл. технологии. 2010. Т. 15, 1. С. 14 32.
77. Агшти А.А., Ильгамоа М.А., Халитоаа Т.Ф. Ударное воздействие струи па жесткую стенку // Актуальные проблемы механики сплошной среды. К 20-летию ИММ КазНЦ РАН. Казань: Фолиант, 2011. Т. 1. С. 134 146.
78. Benjamin Т.В., Ellis А. Т. The collapse of cavitation bubbles and the pressures thereby produced against, solid boundaries // Phil. Trans. R. Soc. Loud. 1966. V. 260, No 1110. P. 221 240.
79. Grum L.A. Surface oscillations and jet. development, in pulsating bubbles // J. de Physique. 1979. V. 40, No 8. P. C8-285 C8-288.
80. Ohl C.-D., Kurz T., Geisler R., Lintlau O., Lauterborn W. Bubble dynamics, shock waves and soiiolumiiiescence // Phil. Trans. R. Soc. Loud. 1999. V. 357, No 1751. P. 269 294.
81. Philipp A., Lauterborn W. Cavitation erosion by single laser-produced bubbles // J. Fluid Mecli. 1998. V. 361. P. 75 116.
82. Аганин A.A., Ильгамоо М.А., Малахов В.Г., Халитоаа Т.Ф., Хисматуллина H.A. Ударное воздействие кавитациоппого пузырька па упругое тело // Учеп. зап. Казап. уп-та. Сер. Физ.-матем. пауки. 2011. Т. 153, кп. 1. С. 131 146.
Поступила в редакцию 26.09.12
Аганин Александр Алексеевич доктор физико-математических паук, профессор. заведующий лабораторией Института механики и машиностроения КазНЦ РАН. E-mail: aganinQkfti.knc.ru
Гусева Татьяна Сергеевна кандидат физико-математических паук, научный сотрудник Института механики и машиностроения КазНЦ РАН. E-mail: gusevaQkfti.knc.ru
