Диссертация (1150462), страница 15
Текст из файла (страница 15)
Для проверки того, что данный эффект не связан с программнойреализацией метода VOF в коде Flag-FS, был также проведен расчет с использованиемпрограммного пакета ANSYS Fluent. Постановка задачи и расчетная сетка были идентичны, дляперевязывания скорости и давления использовался безитерационный метод дробных шагов,шаг по времени был малым (CFL=0,1), для решения уравнения переноса величины Cиспользовался метод геометрической реконструкции (показавший себя в предыдущих тестахлучшим методом из арсенала Fluent). Результаты расчета приведены на рисунке 4.10, б. Видно,что и в данном расчете образовался пристенный слой c заниженными значениями величины C.Если отвлечься от осцилляций поля величины C, имеющих место в данном расчете, можноотметить, что конфигурация слоя сходна с той, что получилась при использовании кодаFlag-FS: тонкий слой вдоль первой половины пути и заполненная воздухом область – вдольвторой.Осцилляциииискривлениеграницслоя,повсейвидимости,связанысдополнительными численными методиками, применяемыми во Fluent и не описанными вдокументации.
Таким образом, можно утверждать, что проблема пристенного воздушного слоя,возникающего при расчете на низкорейнольдсовых сетках, имманентно присуща методу VOF.76В настоящей работе для подавления вышеописанного пристенного воздушного слояприменялась методика дополнительной «положительной» диффузии величины C вблизи стенки(термин «положительная» использован для указания на то, что введенная процедура можетлишь увеличивать значения маркер-функции в ячейках расчетной сетки). Уравнение (2.36)дополнялось членом, отвечающим за диффузию:С C v max C , 0 ,t(4.7)где – коэффициент искусственной диффузии, определяемый по формуле:0d cell d max d cell, d d cell, d d cell.Здесь max – задаваемое пользователем максимальное значение коэффициентаискусственной диффузии, достигаемое на стенке; dcell – расстояние, на котором действуетискусственная диффузия.
Параметр dcell должен быть достаточно малым, чтобы не влиять натечение в целом, но при этом достаточным, чтобы захватить область нефизичного эффекта. Внастоящей работе он задавался равным половине типичного продольного размера ячеек вдольстенки (который примерно равен половинному размеру ячеек в основной части расчетнойобласти).
Параметр max должен быть выше некоторого критического значения, чтобыобеспечить устранение прослойки воздуха (дальнейшее увеличение параметра выше данногозначения уже не влияет на решение). При этом слишком высокие его значения приводят квозрастанию вычислительной сложности решения уравнения (4.7). В настоящей работеиспользовалось значение 1 м2/с.Работает коррекция следующим образом. Если на шаге по времени для ячейки расчетнойсетки вклад от добавленного в уравнение (2.36) диффузионного члена оказываетсяположительным, то он учитывается, в противном случае он не учитывается.
Понятно, чтоподобная схема не является консервативной и приводит к нарастанию общего количествавеличины C (и, следовательно, общего количества жидкости) в расчетной области. Однакофактическиданныйэффектявляетсяпренебрежимомалым:максимальныйобъем«накачиваемой» таким образом жидкости не превосходит (а фактически получается многоменьше) объема пристенного слоя толщиной в dcell.Результаты расчетов с использованием коррекции приведены на рисунках 4.11 и 4.12.Видно, что введенная коррекция обеспечила полное исчезновение нефизичного воздушногослоя, в результате чего на графиках трения о стенку исчезли провалы, и появилось стремлениерешения к сеточно-сошедшемуся при последовательном измельчении сетки.
Также можно77отметить, что в данном течении пристенные функции работают корректно, обеспечивая весьмаблизкие значения трения в широком диапазоне значений y+.Отметим также, что при проведении этих расчетов методика дополнительной диффузииработала одновременно с методикой дополнительного «обострения» межфазной границы,описанной в пункте 3.1.2. При этом данные методики не «мешали» друг другу: методикадиффузии работала только в узкой области вблизи стенки, где она была намного «сильнее»методики «обострения» и быстро устраняла прослойку воздуха.Рисунок 4.11.
Поле величины C вблизи нижней стенки в момент времени 0,4 с: расчет сиспользованием коррекцииРисунок 4.12. Распределение величины трения на нижней стенке: расчеты на сетках сразличными значениями Y+4.1.3. Методика дробных шаговКак было показано в пункте 3.2.2, даже лучшая из исследованных схем аппроксимацииуравнения переноса маркер-функции (2.36) (схема M-CICSAM) требует для обеспеченияприемлемой точности решения проводить расчеты с относительно небольшими шагами повремени, соответствующими числам Куранта до 0,75. Однако при численном решенииуравнений движения жидкости (2.3) и (2.6) можно, как правило, использовать много бóльшие78шаги по времени при сохранении высокой точности решения, экономя тем самымвычислительные ресурсы.
В настоящей работе реализована методика, позволяющая совершатьвнутри одного шага по времени, выполняемого при численном решении уравненийгидродинамики (2.3) и (2.6), несколько более мелких шагов по времени для уравнения переносамаркер-функции (2.36), обеспечивая тем самым оптимальные значения чисел Куранта для всехуравнений. В представленных далее по тексту расчетах один шаг для уравнений (2.3) и (2.6)разбивался на такое количество равных шагов для уравнения (2.36), при котором число Курантана каждом из них не превосходило предустановленного значения 0,5.Отдельное внимание было уделено вопросу, какое поле скорости использовать накаждом из дробных шагов решения уравнения (2.36).
Рассматривались два варианта: (1)линейно интерполировать скорость на текущий «дробный» шаг по значениям с нового истарого «больших» шагов; (2) использовать на всех дробных шагах скорость, полученную какполусумму значений с нового и старого «больших» шагов. Тестирование данных вариантовпоказало, что при одинаковых шагах по времени они обеспечивают примерно одинаковуюточность, поэтому в большинстве расчетов использовался второй вариант как требующийменьшего количества вычислений.4.2.Метод решения уравнений гидродинамики4.2.1.
Вычисление конвективных потоковПричисленномрешенииуравнениябалансаимпульса(2.3)(или2.12),аппроксимированного по методу конечных объемов (2.20), требуется вычисление потоковколичества движения через грани ячеек. Конвективная составляющая этих потоков может бытьсосчитана как f F f v f , где, как и ранее, F f – объемный расход среды через грань f. Посколькуплотность среды в окрестности границы жидкость-газ изменяется на несколько порядков,различные способы аппроксимации ее значений на грани ячеек могут давать значения f ,отличающиеся на порядки, оказывая сильное влияние на получаемое решение.Конечно-объемная аппроксимация уравнения (2.3) предполагает использование такойаппроксимации для значений плотности на гранях ячеек, при которой будет выполнено условиесохранения массы (2.33) (в настоящей работе уравнение (2.33) не решается, а за сохранениемассы в расчетной области отвечает уравнение переноса величины С (2.37), являющеесялинейно-зависимым с ним).
Данное требование будет выполнено при использованиисогласованных аппроксимаций для конвективных потоков величины С в уравнении (2.37) ипотоков количества движения f F f v f , а также для производных по времени в данных79уравнениях, как это делалось, к примеру, в работе [127]: для аппроксимации производной повремени в данных уравнениях использовалась одна и та же схема, а при вычислении потоковимпульса через грани ячеек использовались значения плотности, вычисленные через значениявеличины С на тех же гранях по соотношению (2.4).В настоящей работе по ряду причин было решено отказаться от данной концепции.Одной из таких причин было использование алгоритма, предусматривающего несколькодробных шагов для решения уравнения (2.36) на каждом шаге решения уравненийгидродинамики,описанноговпункте4.1.3.Крометого,использоваласьметодикадополнительного «обострения» фронта, подробно описанная в пункте 4.1.1.Для избавления от упомянутой выше необходимости согласовывать аппроксимацииуравнений, можно переписать уравнение (2.3) в неконсервативной форме (4.8а) (как, например,в работе [135]) или в форме (4.8б) (как в работе [52]).
Форма (4.8б) следует из (2.3) с учетом(2.2) и (2.6). Дискретные аналоги данных формулировок представлены соответственноформулами (4.9а) и (4.9б). В обеих формулировках плотность не входит под знак производной,так что вопрос о способе интерполяции плотности на грань ячейки отпадает.Однако, как выяснилось в ходе методических расчетов, даже в случаях относительнопростых течений использование и той, и другой форм записи приводит к сильной деформациимежфазной границы. В качестве примера, на рисунке 4.13 приведены результаты решениямодельной двумерной задачи о свободном падении круглого «пятна» жидкости в покоящемсявоздухе под воздействием силы тяжести. Расчеты велись с достаточно малым шагом повремени (t = 0,0005 с, число Куранта CFL<0,25) чтобы считать погрешность аппроксимациипо времени пренебрежимо малой.
Тем не менее, для обоих вариантов (4.9), уже послеперемещения пятна на расстояние, меньшее собственного диаметра, обнаруживаются заметныенефизичные искажения его формы, что обусловлено неточным соблюдением баланса импульсав ячейках расчетной сетки. v v v RHStv v v RHStppv ptv pt(4.8а)(4.8б)V p v p F f RHS(4.9а)V p F f v f RHS(4.9б)ff80Здесь RHS – Right-hand Side – обозначает правую часть в уравнениях и включаетисточниковые и диффузионные слагаемые.Рисунок 4.13.
Искажение формы падающего круглого пятна: 1 – расчет по схеме (4.9а), 2 –(4.9б), 3 – (4.10б);фоновой заливкой показано неискаженное пятноВ настоящей работе для ухода от требований, предъявляемых к способу интерполяцииплотности, использовался оригинальный подход, в котором уравнение движения записано вформе (4.10а), получающейся путем вычитания уравнения (2.2), умноженного на скорость, изуравнения (2.3). Дискретный аналог полученного уравнения дается соотношением (4.10б). v v v v v RHStpv ptV f v f F f v p f F f RHSf(4.10а)(4.10б)fБлагодаря такому вычитанию компенсируются ошибки аппроксимации, связанные снеточным выполнением условия баланса массы для ячейки (при этом в аппроксимацияхвторого и третьего слагаемых нужно использовать одинаковые значения плотности на граняхячеек).Результаты расчетов с использованием схемы (4.10б) приведены на рисунке 4.13(пунктирная линия).