Диссертация (1150462), страница 8
Текст из файла (страница 8)
Ее можнорассматривать как обобщенную схему второго порядка.UGD f CD CDf ff ,(2.24)где – задаваемый весовой коэффициент.При = 0,5 схема становится эквивалентной квадратичной интерполяции по тремзначениям A, D и U (см. рисунок 2.2), предложенной в работе [73] и названной QUICK –Quadratic Upstream Interpolation for Convective Kinematics. Схема обладает существеннобольшей устойчивостью, чем схема CD, имеет наименьшую погрешность аппроксимации извсех возможных схем (2.24) для различных значений .При = 2,0 получается схема, эквивалентная линейной экстраполяции по значениям D иU (ячейка U находится вверх по потоку от ячейки D). Cхему далее по тексту будем называтьSOUD – Second Order Upwind Differencing.
Ее погрешность аппроксимации в три раза больше,чем у схемы UGD.Рисунок 2.2. Различные способы аппроксимации величины на грань f по ее значениям вцентрах ячеек U и D, расположенных вверх по потоку и ячейки А, расположенной вниз попотокуСледует отметить, что при использовании противопоточных схем второго порядкааппроксимации возможно развитие нефизичных осцилляций решения, которые могут привестик неприемлемым результатам (например, к появлению отрицательных значений кинетической34энергии турбулентности). Одним из способов борьбы с нежелательными осцилляциямиявляется введение специальных ограничителей, которые в «опасных» местах понижаютпорядок точности схемы, делая ее монотонной.
Широко используется, в частности,ограничитель minmod (2.25), который особенно эффективен в сочетании со схемой SOUD.modmin D R f D ; CDff D ,R u, v uvuv 1 2 uv(2.25)2.2.3. Схемы аппроксимации производной по времениВычисление производной по времени, стоящей в левой части выражения (2.20), требуетиспользования значений дифференцируемой величины с нескольких слоев по времени.Перепишем выражение (2.20) в следующем виде: Pt P flux f n f S ff(2.26)VПраваячастьэтоговыражениясодержитаппроксимациипространственныхпроизводных и не содержит производных по времени. Обозначим ее (). Тогда уравнение(2.20) можно записать в виде pt ()(2.27)Двухслойная схема аппроксимации по времени может быть построена двумя способами: np np1t n np np1t n n1 1 n ,(2.28а) n1 1 n (2.28б)Здесь верхние индексы относятся к слоям по времени (см.
рисунок 2.3); – весовойкоэффициент. Если = 1, то аппроксимация по времени является явной, в противном случаеаппроксимация по времени является неявной. Если = 0,5, то аппроксимация по времени имеетвторой порядок (это схема Кранка-Николсон), в противном случае аппроксимация по времениимеет первый порядок.Трехслойная схема второго порядка с пространственным оператором на последнем слоедля одинаковых шагов по времени имеет вид:3 np 4 np1 np22 tn n(2.29)35Ее обобщение на случай переменных шагов выглядит так: np 2 np1 1 np22 1 tn n(2.30)Здесь индексы n, n-1 и n-2 относятся к слоям во времени (см.
рисунок 2.3), tn 1t n.Рисунок 2.3. Шаблон схемы (2.30) (слева) и схем (2.28) (справа). На шаблоне вертикальный рядточек обозначает время, горизонтальный ряд точек – состояние системы, их пересечениепоказывает, к какому моменту времени относится вычисленное состояние2.2.4. Форма записи уравнения неразрывности и уравнения переноса маркерфункцииИнтегральная форма (2.19), в соответствии с теоремой Остроградского-Гаусса,эквивалентна консервативной дифференциальной форме: flux t(2.31)Консервативную форму имеют уравнения (2.2) и (2.3), в то время как уравнение (2.1)консервативным не является, и не может быть формально использовано в методе конечныхобъемов.
В качестве возможного подхода, можно составить дискретный аналог уравнения (2.1)в следующем виде:C PV vP C f n f S f 0tf(2.32)Однако, использование такой формы может привести к проблемам. Чтобы этопроиллюстрировать предположим, что аппроксимации в уравнениях (2.1) и (2.2) одинаковые, тоесть, что плотность на грани ячеек интерполируется так же, как и величина C.Запишем дискретный аналог уравнения (2.2): PV f v f n f S f 0tf(2.33)Подставляя уравнение (2.4) в уравнение (2.32), получим дискретное соотношение (2.34),в котором используются те же аппроксимации, что и в выражении (2.33):36 PV vP f n f S f 0tf(2.34)Вычитая уравнение (2.34) из уравнения (2.33), получим, что численно будет выполненоследущее соотношение: vfP v f n f S f 0(2.35)fАнализируя это соотношение, можно прийти к выводу, что в поле скорости будутвозникать сильные нефизичные флуктуации вблизи межфазной границы, где имеет местоизменение плотности на несколько порядков.
Если на разных гранях ячейки плотностьразличается во много раз, то небольшое отличие между скоростями в центре ячейки и на граняхс большой плотностью потребует огромных отличий между скоростями в центре ячейки и награнях с малой плотностью, что означает большие нефизичные различия в скоростях соседнихячеек.Болееоптимальнымподходомпредставляетсяприведениеуравнения(2.1)кконсервативной форме, что может быть сделано при несжимаемости жидкости и газа (2.6):С C v 0t(2.36)При этом в силу линейной связи (2.4) между плотностью среды и маркер-функций Cуравнение (2.36) является линейно-зависимыми с уравнением неразрывности (2.2). Чтобыопределяющие уравнения были линейно-независимыми, вместо уравнения (2.2) используетсяусловие несжимаемости среды (2.6).Для уравнений (2.36) и (2.6) дискретные конечно-объемные аналоги имеют следующийвид:С ptFV C f Ff 0(2.37)ff0(2.38)fЗдесь, как и ранее, F f S f n f v f – объемный расход через грань f.Решение уравнения для величины C в форме (2.36) обеспечивает баланс величины C (и,как следствие, баланс массы) в расчетной области в целом.373.
Тестирование известных численных схем для решения уравненияпереноса маркер-функции3.1.Описание схем аппроксимации уравнения переноса маркер-функции3.1.1. Предварительные замечанияКак было сказано ранее, в методе VOF для определения положения свободнойповерхности необходимо решать уравнение (2.36) конвективного переноса объемной долижидкости. Корректное решение данного уравнения является непростой задачей, посколькунеобходимо обеспечить одновременное выполнение двух условий: удерживать значениямаркер-функции внутри интервала 0-1 и не допускать «размытия» контактной границы (то естьобласти, где имеет место изменение маркер-функции от 0 до 1) на множество ячеек расчетнойсетки.
Выбор подходящей схемы аппроксимации уравнения (2.36) предполагает детальноесравнение возможностей различных схем.Среди пользующихся популярностью в настоящее время численных схем для решенияуравнения переноса маркер-функции можно отметить т.н. «сжимающие» схемы CICSAM [126]и HRIC [91]: они реализованы во многих коммерческих CFD-кодах. В программном пакетеANSYS Fluent реализован также т.н. метод геометрической реконструкции, которыйрекомендуется для решения уравнения (2.36). Кроме того, относительно недавно появиласьпривлекательная схема M-CICSAM [133].Необходимостьвразработкеспециализированных«сжимающих»схембылаобусловлена тем, что «стандартные» противопоточные схемы, традиционно применяемые дляаппроксимации конвективных слагаемых, как правило, не обеспечивают необходимогокачества решения уравнения переноса маркер-функции.Тестовые расчеты показывают, что «стандартные» схемы вида (2.24) приводят квозникновению «вылетов» значений маркер-функции C за границы допустимого интервала 0–1.Одна из первых попыток построить схему, одновременно удерживающую межфазную границуот «размытия» и значения маркер-функции в интервале 0–1 была предпринята в работе [51] приразработке метода VOF.
Авторы работы рассматривали схемы первого порядка саппроксимацией вверх и вниз по потоку (далее по тексту UD и DD), то есть сводящиеся кприравниванию значения на грани f к значению соответственно в ячейке D или A (см.рисунок 2.2). Известно, что данные схемы обладают сильной численной диффузией, причемсхеме DD соответствует диффузия с отрицательным коэффициентом («антидиффузия»).Применительно к решению уравнения (2.36) это означает, что схема UD приведет к сильному38«размытию» межфазной границы, а схема DD благодаря обострению градиентов сделаетграницу «острой», но при этом, приводя к увеличению локальных максимумов и уменьшениюлокальных минимумов, не обеспечит ограниченности решения. В работе [51] предпринятапопытка устранить недостатки схем DD и UD путем поиска оптимального их сочетания.Значения маркер-функции С на гранях ячеек предлагалось определять как линейнуюкомбинацию ее значений из соседних ячеек D и A, причем веса значений в данной линейнойкомбинации определялись исходя из самих значений в ячейках D и A.
Предложенная схема необеспечивала полного отсутствия «вылетов», а лишь ограничивала их до относительнонебольших значений. Авторы предлагали «обрезать» зашкалившие значения (принудительно накаждом шаге по времени приравнивать значения величины C к единице в ячейках, где ониоказались выше единицы, и к нулю, где ниже нуля). Таким образом, оригинальная схема дляметода VOF [51] формально не являлась консервативной и не обеспечивала сохранение массыжидкости в объеме. Кроме того, схема приводила к искусственному искривлению межфазнойграницы, когда вдоль нее направлена скорость потока [51, 69, 74] (в частности, плавные волны,бегущие по поверхности жидкости, обострялись до «ступенек»).
Для подавления этого эффектав работе [51] предложено переключаться на полностью противопоточную схему UD, когда уголмежду межфазной границей и направлением потока меньше 45о. Однако такой подход оказалсяслишком грубым: результаты расчетов существенно зависели от предустановленной величиныугла, на котором производилось переключение на противопоточную схему [69], и, кроме того,не удавалось полностью избавиться от искривления межфазной границы [11, 69, 108, 127].3.1.2. Диаграмма нормализованной переменной (NVD) и критерий локальнойограниченности (CBC)В работе [35] для упрощенного случая одномерного течения формально показано, чтодля обеспечения ограниченности решения необходимо учитывать значение величины C вячейке, находящейся вверх по потоку от ячейки D (ячейка U, см. рисунок 2.2).