ChISLENNOE MODELIROVANIE NESTATsIONARNYKh TEChENIY NESZhIMAEMOY ZhIDKOSTI SO SVOBODNOY POVERKhNOST'Yu PRI POMOShchI VOF METODA

Abstract

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

Full Text

Течения несжимаемой жидкости со свободной по- верхностью очень часто встречаются как в различных природных явлениях, так и во множестве технологи- ческих процессов. Примеры таких течений можно найти практически в любых промышленных объектах. Поэтому умение эффективно решать подобного рода задачи является очень важным в процессе проектирова- ния и оптимизации различных устройств и аппаратов. На сегодняшний день существует множество алго- ритмов решения задач со свободной границей. Это, прежде всего, так называемые body-fitted методы [1], в которых свободная поверхность жидкости отслежива- ется узлами расчетной сетки. Естественно, что при таком подходе расчет перемещения жидкости в про- странстве требует пересчета расчетной сетки на каж- дом временном шаге, что может быть весьма затрат- ным. Кроме того, поскольку форма свободной по- верхности жидкости и ее движение часто очень слож- ны, то использование body-fitted методов может при- вести к существенному искривлению расчетных ячеек, что приводит к дополнительной погрешности в ре- зультатах. В последнее время, в связи с возросшей вычисли- тельной мощностью, для решения задач со свободной поверхностью стали активно использоваться так на- зываемые бессеточные методы или методы частиц [2]. Основная идея этих методов состоит в том, что жид- кость представляется конечным числом взаимодейст- вующих частиц, движение которых подчиняется урав- нениям движения. В этих методах не требуется раз- биение расчетной области на конечное число ячеек. Хотя их идеология продолжает интенсивно развивать- ся, их реальное применение для решения сложных пространственных задач из-за больших вычислитель- ных затрат пока весьма ограничено. Обзор литературы показывает, что хотя на сегодняшний день существует огромное количество под- ходов к решению задач со свободной границей, самым эффективным и распространенным по применению является VOF метод (Volume of Fluid) [3] (и его более современные модификации). Данный метод сочетает в себе высокую вычислительную эффективность и дос- таточную точность разрешения межфазной границы. Благодаря своей простоте он относительно легко реа- лизуется в рамках алгоритма расчета несжимаемых течений жидкости, основанного на решении системы уравнений Навье-Стокса. К достоинствам метода также стоит отнести возможность описывать такое важное для течений со свободной поверхностью явле- ние, как поверхностное натяжение. Моделирование данного явления достаточно просто проводить при помощи VOF метода с использованием CSF модели (continuum surface force) [4]. Учитывая все преимущества VOF метода, авторы работы реализовали его в программном комплексе «sFlow» в качестве алгоритма решения задач со сво- бодной поверхностью. Математическая модель движения несжимае- мой жидкости со свободной поверхностью. Течение двухфазной среды в ламинарном случае описывается системой уравнений Навье-Стокса. Ниже приведены уравнения данной модели, выражающие основные законы сохранения. Уравнение неразрывности (закон сохранения масв введении в уравнение количества движения (2) объ- емной силы FS - величина которой находится из соот- ношения FS = skÑF , где k - кривизна свободной поверхности, определяе- мая как дивергенция вектора нормали: æ n ö image ç n ÷ k = Ñ ç . è ø Нормаль к свободной поверхности вычисляется, в свою очередь, как градиент объемной доли жидкой фазы в ячейке: n = ÑF. На твердой стенке величина вектора нормали оп- ределяется по краевому углу смачивания q : сы двухфазной среды): n = n cos(q) + t sin(q), image ¶r + Ñ(rv) = 0. ¶t (1) w w где nw, фw - нормальный и тангенциальный к стенке вектора. Уравнение количества движения (закон сохране- ния импульса): ¶rv Численный алгоритм и программная реализа- ция модели. Для компьютерной реализации матема- тической модели использовался разрабатываемый авторами программный комплекс «sFlow». Пакет проimage ¶t + Ñ(rv × v) = -Ñp + Ñ(t) + rg + FS , (2) грамм «sFlow» позволяет моделировать стационарные где FS - сила поверхностного натяжения; g - ускоре- ние свободного падения; ф - тензор вязких напряже- ний, составляющие которого определяются по выра- жению и нестационарные, ламинарные и турбулентные тече- ния несжимаемых газов, одно и многокомпонентных неизотермических запыленных сред с химическими реакциями и сложным теплообменом, со свободной поверхностью. При расчетных исследованиях различ- æ ¶u ¶u j 2 t = m i + - d ¶uk ö. ных объектов возможна настройка пакета программ на image image image ij ç ¶x ¶x 3 image ij ¶x ÷ расчет соответствующих процессов. Пакет программ è j i k ø Для определения положения свободной поверхно- сти жидкости в пространстве, как уже было сказано, будем использовать VOF метод [3]. Идея метода со- стоит в том, что поскольку свободную поверхность можно рассматривать как границу раздела двух фаз (газ-жидкость), то каждую фазу можно задать своей объемной долей. Пусть F(x, y, z, t) - доля жидкой фазы в расчетной ячейке, которая определяется следующим образом: ì0, если ячейка пустая; F (x, y, z, t) = í î1, если ячейка заполнена жидкостью, и 0 < F(x, y, z) < 1, если через ячейку проходит граница раздела фаз. Отслеживание перемещения свободной границы осуществляется путем решения уравнения переноса для объемной доли: включает в себя модуль построения геометрии и рас- четных сеток, систему визуального анализа простран- ственных полей и постпроцессинга, модуль решения уравнений гидродинамики. Разностный аналог конвективно-диффузионных уравнений находится с помощью метода конечного объема для структурированных многоблочных сеток [5; 6], при применении которого автоматически вы- полняется консервативность полученной схемы. Для аппроксимации конвективных членов уравнений гид- родинамики (1), (2) используется противопоточная схема второго порядка QUICK [5]. Для аппроксима- ции нестационарных слагаемых уравнений гидроди- намики используется неявная схема первого порядка. Диффузионные потоки и источниковые члены ап- проксимируются со вторым порядком точности. Связь между полями скорости и давления реализуется при помощи SIMPLEC процедуры на совмещенных сетках. Для устранения осцилляций поля давления исimage ¶F + v × ÑF = 0. ¶t (3) пользуется подход Рхи-Чоу, связанный с введением монотонизатора в уравнения для поправки давления Плотность и молекулярная вязкость двухфазной среды находятся через объемную долю жидкости в ячейке по правилу смеси: r = r1 F + (1- F )r2 , m = m1 F + (1- F )m2 . Для моделирования поверхностного натяжения применялся CSF алгоритм [4]. Суть алгоритма состоит [5; 7]. Полученная система разностных уравнений ре- шается итерационным способом с применением мето- дов неполной факторизации (явный метод Булева, ускоренный методом сопряженных градиентов) [8]. Численные эксперименты показали, что на достоверность результатов, полученных при использовании VOF метода, существенное влияние оказывает качест- во метода решения уравнения (3). Для интегрирования уравнения (3) авторами были рассмотрены следующие методы: явная TVD схема Superbee [5], явная TVD схема Superbee [5] с локально одномерным расщепле- нием по пространству [9], неявная схема TVD схема Superbee первого и второго порядка аппроксимации по времени, схема первого порядка UDS [6], противо- поточная схема второго порядка QUICK. Верифика- ция данных схем аппроксимации проводилась на множестве как одномерных, так и пространственных конвективных задач. По итогам тестовых расчетов лучшей для решения уравнения (3) оказалась явная TVD схема Superbee c использованием локально од- номерного расщепления. Все представленные ниже результаты расчетов получены при помощи данной схемы. Тестирование расчетного алгоритма. Задача об обрушении водяного столба. Данная задача является каноническим тестом для алгоритмов решения задач со свободной поверхностью. Постановка и начальное условие этой двумерной задачи представлены на рис. 1. Стенка высотой 2а ограничивает столб воды шири- ной а, в начальный момент расчета стенка мгновенно удаляется и под действием силы тяжести жидкость растекается по расчетной области. Для расчета ис- пользовалась прямоугольная однородная сетка, со- стоящая из 100 Ч 100 узлов. Число Куранта задавалось равным 0,5. Наглядное представление о положении свободной поверхности воды в пространстве в раз- личные моменты времени дает рис. 2 (черный цвет - вода). Ситуация, когда вода дошла до правого края расчетной области и стала натекать на стенку, пред- ставлена на рис. 2, г. Данные эксперимента на этот случай уже не распространяются. Количественное сравнение с экспериментом [10] представлено на рис. 3. Сопоставление с эксперимен- том проводилось по двум параметрам: расстоянию x (рис. 1), на которое распространится вода за время t и высоте водяного столба b за это же время. Наблюдаетimage ся достаточно хорошее совпадение численных резуль- татов с экспериментальными данными (рис. 3) [10]. 2a b а х Рис. 1. Схема задачи о разрушении дамбы Моделирование процесса заливки жидкого ме- талла в изложницу [11]. Приведем пример использо- вания реализованного для оптимизации заливки жид- кого металла в формы в процессе изготовления алю- миниевых слитков. Важным условием повышения эффективности производства алюминиевых слитков является условие увеличения скорости разлива метал- ла в изложницы. Но при этом неизбежно происходит увеличение потерь металла, обусловленное разбрыз- гиванием и окислением кислородом воздуха. В связи с этим было проведено исследование влияния скорости разливки и способа литья на качество получаемых слитков. Было рассмотрено два различных способа разливки жидкого алюминия: разливка при помощи разливочной машины (струя металла падает верти- кально) и разливка при помощи лотка (струя падает под углом к дну изложницы). Расчеты проводились для трех различных скоростей литья: 5, 7 и 10 т/ч. t = 0,15 c t = 0,25 c a б t = 0,3 c t = 0,4 c image в г Рис. 2. Динамика разрушения дамбы Эволюция формы свободной поверхности жидкого металла представлена на рис. 4 (разливка с лотка, ско- рость литья - 7 т/ч). Видно, что процесс литья сопро- вождается значительным волнообразованием. Более того, в некоторых вариантах при высоких скоростях литья (более 12 т/ч) наблюдалось выплескивание ме- талла за пределы изложницы, что является недопус- тимым. Кроме того, было обнаружено такое интерес- ное явление, как захват пузырей воздуха падающей струей металла. Наличие дополнительных пузырей воздуха под свободной поверхностью металла приво- дит к увеличению площади контакта алюминия с ки- слородом, а следовательно, и к ухудшению свойств слитков. В качестве относительной «количественной» ха- рактеристики качества разливки, использовался такой критерий, как избыточная площадь поверхности жид- кого металла. Под избыточной площадью поверхности здесь понимается разность между «реальной» площа- дью поверхности металла в данный момент времени и площадью поверхности металла в случае идеальной разливки (без волнообразования и попадания пузырей воздуха). Чем выше величина избыточной площади поверхности металла, тем выше площадь соприкосно- вения металла с кислородом воздуха (рис. 5). С увели- чением скорости литья происходит увеличение избыточной площади поверхности жидкого металла. Но поскольку при повышении скорости литья пропор- ционально уменьшается и время разливки (а значит и время контакта жидкого металла с воздухом), то мож- но сказать, что существует некая оптимальная вели- чина скорости литья, при которой среднее по времени значение избыточной площади жидкого металла будет минимальным (а следовательно, минимальны будут относительные потери металла на окисление). В дан- ном варианте разливки оптимальная скорость литья близка к 10 т/ч. По итогам проведенного исследования обнаружено, что разливка с помощью разливочной машины приводит к более интенсивному процессу волнообра- зования и захвату пузырей газа струей жидкого ме- талла, что влечет увеличение потерь алюминия на окисление. Поэтому вариант разливки алюминия при помощи лотка является более предпочтительным. Увеличение скорости заливки металла с 5 до 10 т/ч, хотя в целом и приводит к увеличению избыточной поверхности алюминия, дает уменьшение средней по времени площади контакта кислорода воздуха с жид- ким металлом. Данная величина может служить каче- ственным показателем величины потерь алюминия на окисление. В связи с этим можно рекомендовать уве- личение скорости литья с 5 до 10 т/ч. t g / a t g / a image image Рис. 3. Зависимость безразмерных размеров столба от безразмерного времени а б в Рис. 4. Распределение свободной поверхности металла в различные моменты времени: а - t = 2 c; б - t = 7 c; в - t = 9 c 0,14 0,12 0,1 0,08 0,06 0,04 0,02 10 т/ч 5 т/ч 7 т/ч t, c image Рис. 5. Зависимость избыточной поверхности жидкого металла от времени заливки и скорости литья Таким образом, реализован эффективный алгоритм расчета нестационарных течений несжимаемой жид- кости со свободной поверхностью с учетом сил по- верхностного натяжения на основе VOF метода. Про- ведено решение ряда тестовых задач, в ходе чего по- лучено хорошее качественное и количественное сов- падение численных результатов с экспериментальны- ми данными. Реализованный алгоритм был успешно применен к решению задачи оптимизации разливки жидкого металла по изложницам. Полученные резуль- таты позволили сделать предложения, которые спо- собствовали существенному повышению эффективно- сти исследованного технологического процесса.