Диссертация (1150462), страница 16

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 16)

Видно, что искажения формы пятна практически отсутствуют. При этомодной линией отображены практически идентичные результаты двух расчетов с разнымиспособами вычисления значений плотности на гранях ячеек: по соответствующим значенияммаркер функции C (как это необходимо при решении уравнения в форме (2.3)) и как среднеегеометрическое плотностей в центрах смежных ячеек. Таким образом, дискретизация,основанная на записи уравнения сохранения импульса в форме (4.10а), обеспечивает высокоекачество решения и его малую чувствительность к способу вычисления плотности на граняхячеек.81Предлагаемая форма записи для уравнения движения естественным образом переноситсяна уравнение сохранения произвольной скалярной величины, которое можно представить вформе обобщенного уравнения конвективно-диффузионного переноса:   v      v   RHSt(4.11)Здесь  – переносимая величина (например, компонента скорости ui или параметртурбулентности), RHS – правая часть уравнения, которая может включать источниковые члены,диффузионные слагаемые и др.Дискретный аналог уравнения (4.11) имеет вид:p ptV    f  f F f   p   f F f  RHSf(4.12)fОтметим, что в настоящей работе для вычисления значений величины  на гранях ячеекиспользовались схемы UGD и SOUD (первая – для компонент скорости, вторая – дляпараметров турбулентности) с дополнительным ограничителем minmod (2.25).4.2.2.

Вычисление диффузионных потоковДиффузионный перенос имеет место в уравнениях переноса параметров турбулентности(2.9), (2.10), (2.14), (2.15), а также в уравнении движения (2.3) (или 2.12). В первых четырехуравнениях за эффекты диффузии отвечают слагаемые вида   ,(4.13)где  – коэффициент диффузии (в общем случае неоднородный в пространстве),  –переносимая величина (k,  или ).В уравнении движения за эффекты диффузии (вязкости) отвечает слагаемоет  τ     v  v  .(4.14)Здесь  обозначает турбулентную вязкость.Дискретные аналоги данных членов при использовании метода конечных объемовимеют соответственно вид:Sf f n f   ff(4.15)VSfτ f nffVSf f n f  v  ffVSfт f n f  v  ffV(4.16)82Выражение (4.15) и первое слагаемое в выражении (4.16) могут быть преобразованы квидуSff    f nf   n  fV(4.17)Здесь в роли  также могут дополнительно выступать и компоненты скорости.Значения коэффициентов диффузии  на гранях ячеек вычисляются по схеме (2.21).

Под  величиной   здесь понимается производная от величины  вдоль нормали к грани ячейки. n  fВ коде Flag-FS для вычисления значения этой производной с учетом возможной скошенностиячеек используется следующее преобразование:          n f   f   f , n  f    f(4.18)где  f – единичный вектор, направленный от центра P рассматриваемой ячейки к центру  ячейки M, соседствующей с ней по грани f (см. рисунок 4.14),   – производная от   fвеличины  по этому направлению.PfMnРисунок 4.14.

К расчету диффузионного потока через грань ячейки    PВеличина   вычисляется по формуле M. Величина градиента от величины    fPMна грани  f вычисляется по схеме (2.21) как линейная интерполяция значений градиента вцентрах ячеек P и M.При вычислении значения второго слагаемого в выражении (4.16) значения компонентттензора v  f на грани ячейки вычислялась путем линейной интерполяции значений этихкомпонент из центров ячеек P и M, вычисляемых с использованием выражения (2.23).

Отметим,83что в условиях несжимаемой среды и с учетом того, что в большей части расчетной области (заисключением областей вблизи межфазной границы) эффективная вязкость меняется плавно,вклад от второго слагаемого в выражении (4.16) является сравнительно малым.4.2.3. Аппроксимация по времениДля выбора схемы аппроксимации по времени применительно к уравнению (4.11)сравнивались две схемы аппроксимации второго порядка: неявная трехслойная схема (4.19) исхема Кранка-Николсон (4.20): n  2      n 1 1      n  2    n  n v n   n    n v n  RHS  n 1  2n n   n 1~ v~  ~    v~  RHS ~nt (4.19)(4.20)Здесь верхние индексы относятся к слоям во времени (см. рисунок 4.15),   t n -1 t n .Рисунок 4.15.

Шаблон схемы (4.19) (слева) и схемы (4.20) (справа)~Величины  , ~ и v в схеме (4.20) относятся к промежуточному слою tn-1/2. Значениеплотности на этом слое вычисляется как полусумма значений на слоях n-1 и n, а в выражение~ (и аналогично для v~ ) введена оригинальная стабилизирующая поправка,для предотвращающая развитие осцилляций во времени: n 1   n,2nn 1~      K    n2 1   n1   n 1 21    1   (4.21)Как можно видеть, введенная поправка пропорциональна второй производной  2  t 2 иимеет второй порядок малости.

Коэффициент K должен быть достаточным для обеспеченияустойчивости численной схемы, но не слишком большим, чтобы не искажать получаемоерешение. В настоящих расчетах принималось значение K=0,125.В качестве теста для сравнения схем (4.19) и (4.20) рассматривалась уже знакомая (см.пункт 4.2.1) тестовая задача о падении «пятна» в форме круга в поле силы тяжести. Расчетыпроводились с переменным шагом по времени, обеспечивающим постоянное значение числа84Куранта, которое в различных расчетах задавалось равным 0,2, 0,4 и 0,8 (при CFL=0,8 накаждый шаг решения уравнений гидродинамики приходились два дробных шага для уравнения(2.36)).

Длительность падения составляла 0,12 с, в течение которых пятно успевалопереместиться на расстояние более 3,5 своих диаметров (см. рисунок 4.16). Конечная формапятна приведена на рисунке 4.17.симметрияРисунок 4.16. Постановка задачиРисунок 4.17. Конечная форма падающего круглого пятна, рассчитанная с помощью схемы(4.20) (слева) и (4.19) (справа); фоновой заливкой показано неискаженное пятноВидно, что схема (4.20) способна сохранять форму пятна практически неискаженной причислах Куранта до CFL=0,4, а при CFL=0,8 искажения не превосходят половины размераячейки сетки (за исключением небольшой области в окрестности лобовой точки).

Прииспользовании схемы (4.19) искажения заметны уже при CFL=0,4, а при CFL=0,8 искаженияуже весьма существенны. В целом, можно сказать, что схема (4.20) обеспечивает приемлемуюточность при числах Куранта до 0,8, а схема (4.19) требует почти в 2 раза более мелкого шагапо времени. Также в ходе тестов обнаружено, что схема (4.19) требует больше итераций накаждом шаге по времени. Таким образом, схема (4.20) существенно экономит вычислительныересурсы, в связи с чем она была выбрана для использования в последующих расчетах.4.2.4.

Вычисление градиента давленияПри использовании метода конечных объемов аппроксимация градиента давления вуравнении баланса импульса определяется выражением:p P  1  p fVfS f nf(4.22)85Обычно для вычисления значений pfна гранях ячеек используют линейнуюинтерполяцию (4.23) (представленные выше результаты получены с ее помощью), однако онаработает некорректно, если на границе жидкость-газ имеет место разрыв градиента давления.p f lin p1d 2  p2 d1d1  d 2(4.23)Пояснить это можно на примере, в котором граница жидкости, покоящейся в поле силытяжести, совпадает с границей между ячейками сетки (рисунок 4.18). В этом случае приправильных значениях давлений p1 и p2 в центрах примыкающих к границе ячеек, значение plinна грани f будет завышенным. Как следствие, в численных аналогах уравнения балансаимпульса (2.3), записанных для ячеек 1 и 2, аппроксимации для градиента давления и объемныхсил оказываются не согласованными: сумма сил давления по всем граням не будет равнадействующей на ячейку силе тяжести.

Это не позволит получить корректное решение(например, в окрестности фронта могут развиваться осцилляции скорости и давления, вплоть доразвала численного решения).Рисунок 4.18. Распределение давления в покоящейся системе жидкость-газВ программном пакете Fluent в качестве одной из опций для интерполяции давленияпредложена схема body-force weighted interpolation, рекомендованная в описании к Fluent длярасчетов течений жидкости со свободной поверхностью. Однако достаточно полное описаниеданной схемы найти в литературе не удалось.Проблема согласования аппроксимаций различных членов в уравнении балансаимпульса при расчете течений со свободной поверхностью по-прежнему остается актуальной, ипоиск ее решения продолжается.

В частности, в недавней работе [81] предложен следующийспособ. Объемная сила, действующая на конечный объем, вычисляется путем осреднениязначений данной силы на гранях контрольного объема, для градиента давления используетсястандартная аппроксимация (4.23). Такой подход гарантированно обеспечивает корректноерешение в случае покоящихся жидкости и газа, когда объемные силы уравновешиваются86градиентом давления. Однако, вопрос о применимости данного подхода в условиях различныхвидов движения жидкости и газа, в частности, при свободном падении капли, когда градиентыдавления малы, а объемные силы уравновешиваются конвективным и локальным ускорением,представляется нетривиальным и требует отдельного систематического исследования, котороене было в полном объеме проведено авторами упомянутой работы. Причем при решении задачио всплытии пузыря в жидкости в работе [81] получено решение с существеннымиосцилляциями скорости вблизи межфазной границы (даже с использованием предложеннойкоррекции).В настоящей работе опробован иной прием избавления от осцилляций скорости вблизимежфазной границы.

Поскольку в задаче с покоящейся жидкостью отношение величинградиентадавлениясразныхсторонотмежфазнойграницыравноотношениюсоответствующих плотностей среды, для получения правильных значений на гранях ячеекможно использовать интерполяцию (4.24) с весами, равными значениям плотности в соседнихячейках. При этом расчетный градиент давления будет в точности компенсировать действиеобъемных сил (см. рисунок 4.18).p f weighted p1 2 d 2  p2 1d11d1   2 d 2(4.24)Рисунок 4.19. Поле давления в окрестности конечного положения падающего круглого пятна(на момент времени 0,12 с), полученное с использованием схемы (4.24) (слева) и схемы (4.25)(справа).

Результаты расчетов на исходной (сверху) и измельченной (снизу) сетках87Однако выяснилось, что при решении рассматриваемой выше задачи со свободнымпадением жидкости схема (4.24) приводит к возникновению сильных четно-нечетныхосцилляций давления вблизи межфазной границы (см. рисунок 4.19, слева-сверху) и кискажению ее формы. Причем данный эффект сохраняется при дроблении расчетной сетки (см.рисунок 4.19, слева-снизу). Проблему удалось решить путем комбинирования схем (4.23) и(4.24) в зависимости от угла  между нормалью к грани ячейки f и силой тяжести:p f    p f weighted  1    p f lin ,  cos 2  (4.25)Схема (4.25) корректно работает как в случае покоящейся жидкости, так и в случаепадения «пятен» жидкости в поле силы тяжести (см. рисунок 4.19, справа).4.2.5.

Характеристики

Список файлов диссертации