Минимальный базис задач валидации методов расчета течений со свободной поверхностью Текст научной статьи по специальности «Физика»

УДК 532.54

А.С. Козелков1'2, А. А. Куркин2, И. Л. Шарипова1, В.В. Курулин1, Е.Н. Пелиновский2, Е.С. Тятюшкина1, Д.П. Мелешкина1, С.В. Лашкин1, Н.В. Тарасова1

МИНИМАЛЬНЫЙ БАЗИС ЗАДАЧ ВАЛИДАЦИИ МЕТОДОВ РАСЧЕТА ТЕЧЕНИЙ

СО СВОБОДНОЙ ПОВЕРХНОСТЬЮ

ФГУП «РФЯЦ-ВНИИЭФ» 1,

2

Нижегородский государственный технический университет им. Р.Е. Алексеева

Цель: В работе обсуждается процесс валидации, являющийся важным этапом на пути промышленного внедрения инженерных пакетов программ.

Метод: В основе исследования лежат физико-математические и численные модели, предназначенные для моделирования течений вязкой несжимаемой жидкости со свободной поверхностью.

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

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

Ключевые слова: валидация, инженерные пакеты программ, течения, свободная поверхность, численное моделирование.

Введение

Данная работа является продолжением работы [1], в которой представлен минимальный базис задач, необходимый для валидации программ (или модулей в составе пакета программ), моделирующих течения вязкой несжимаемой жидкости со свободной поверхностью. Процесс валидации, определяющий степень точности, с которой модель отображает объективную реальность с точки зрения предполагаемого использования, является важным этапом на пути научно-промышленного внедрения разрабатываемых программ. Процессу валида-ции, сопутствует калибровка, в ходе которой осуществляется регулирование числовых или физических параметров модели с целью лучшего соответствия экспериментальным данным.

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

В процессе валидации проверяются научные основы модели путем сравнения с экспериментальными данными и устанавливается, согласуются ли результаты численного моделирования с физическими реалиями. Степень точности, требуемая от результатов моделирования, зависит от предполагаемого их использования. Как правило, к каждой расчетной величине, будь то интегральная или абсолютная характеристика, предъявляются свои требования. Изложение основных принципов верификации и валидации методов CFD-моделирования представлены в [1, 2]. Процедура валидации, состоящая из стадии оценки сходимости итераций, стадии проверки решения на непротиворечивость, стадии сравнения с экспериментом, а также содержащая уровни классов задач представлена в [1].

© Козелков А.С., Куркин А.А., Шарипова И.Л., Курулин В.В., Пелиновский Е.Н., Тятюшкина Е.С., Мелешкина Д.П., Лашкин С.В., Тарасова Н.В., 2015.

Одним из наиболее сложных классов задач вычислительной гидродинамики является моделирование течений с поверхностями раздела двух или более сред, которые в свою очередь можно подразделить на два класса - смешивающиеся и несмешивающиеся жидкости. Задачи этого класса встречаются во многих инженерных и научных задачах и играют большую роль в природе и технике. К этим задачам относят струйные течения, движение пузырей, образование пленок, движение капельных сред и т.д. В общем случае физико-математические модели для моделирования этого класса течений подразумевает использование многофазного подхода, учитывающего межфазное взаимодействие. Такие модели достаточно сложны и должны учитывать специфику моделирования конкретного течения. Отметим здесь, что экспериментальные и натурные исследования таких течений сопряжены со значительными трудностями и затратами, поэтому оценка степени достоверности численных методов, моделирующих такие течения, является очень актуальной задачей.

Распространенным частным случаем таких течений являются течения со свободной поверхностью. К таким течениям относятся течения топлива в форсунках двигателей, движение смазочных материалов, поверхностные волны на воде и др. Для такого класса задач ввиду того, что плотность фаз, например, жидкости и воздуха, на порядки различаются друг от друга, предполагают, что на движение жидкости газ не оказывает особого влияния и ограничиваются предположением о том, что взаимодействующие фазы не смешиваются. Моделирование таких течений осуществляется с помощью специализированных моделей, основанных на выделении объемной доли вещества в объеме [3].

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

1. Метод расчета задач со свободной поверхностью

Будем рассматривать несжимаемые течения, состоящие в общем случае из произвольного количества гетерогенных фаз разделенных границей раздела. В рамках односкоростной модели смеси движение данной системы описывается консервативными уравнениями, включающими уравнения неразрывности, сохранения импульса и уравнение для объемных долей фаз, которые имеют вид [4, 5]:

V • и = 0,

12 а(к V* >и = ■-V • 2 [*(к ¥к }ии)+ V • 2 [*(к к ц(к ^и)-^ + pg, (1)

^ к к к

- а(к }р(к} + V • (а(к }р(к} и )= 0,

здесь и - вектор общей скорости совокупности фаз к, р(к) - плотность фазы к, а(к) - объемная

доля фазы к и 2 а к = 1, Р - давление, д(к) - молекулярная вязкость фазы к, g - ускорение

к

свободного падения.

Основная сложность при численном решении системы (1) заключается в определении связи поля давления с полем скорости. Процедура согласования поля давления с полем скорости должна приводить к одновременному удовлетворению уравнений неразрывности и со-

хранения импульса. Наиболее распространенными являются методы типа SIMPLE, основанные на процедуре коррекции давления или принципе расщепления неизвестных [6, 7].

Для получения SIMPLE-процедуры численного решения системы (1), опустим уравнение для объемных долей и массовые силы, рассмотрим ячейку Р с гранями f = nb(P) и запишем систему уравнений (1) в дискретном виде:

2

u }Sf = О,

f=nb( P)

г

2 a(kV(k)u

n-1

u

-V = - 2 2 a(k )P f )u ^u fSf +

(2)

f=nb(P) k

+ 2 2a Vf )VufSf - 2pfSf +2a )P(k)gV,

f=nb(P) k

f=nb(P)

где п - временной слой, т - временной шаг, - площадь грани/ разделяющей контрольные объемы расчетной сетки (рисунок 1), и/ - величина скорости на грани (здесь и далее индекс / означает принадлежность величины к грани), пЬ(Р) - количество граней ячейки (в данном случае ячейки Р), см. рис. 1.

Наиболее удобно использовать совмещенный алгоритм [8] решения системы (2), для реализации которого систему (2) перепишем в виде:

2

f=nb( P)

n . т\

u f + D f

Vpnf- -Vpnf

Sf=0,

2 a(k )p(k )-

(un - u n-1)

V = -2 2 a(k )pfk )u n-1u }Sf +

f=nb(P) k

(3)

+ 2 2a(kV(k)VufSf - 2PfSf +2a(k)P(k)gV.

f=nb(P) k

f=nb( P)

Здесь «верхняя черта» означает, что величина на грани получена интерполяцией из соседних контрольных объемов. В уравнение неразрывности в системе (3) была использована поправка Рхи-Чоу [9], нивелирующая разницу приближений градиента давления в уравнениях неразрывности и сохранения импульса, а также связывающая поля скорости и давления при одновременном решении уравнений неразрывности и движения. В алгебраической форме система имеет вид:

2 aP jp +

pj p

je{ p,u,v,w}

2

N=NB(p P)

2 aN j N

p

N

je{ p,u,v,w}

2aPjP + 2

je{p,u,v,w} N=NB( P)

= bp,

= bp, i = {u, v, w}.

(4)

X

Суммирование по индексу «у» для первого уравнения системы (4) - уравнения неразрывности, дает коэффициенты общей матрицы системы для вычисления давления в контрольных объемах дискретной модели. Данные коэффициенты имеют вид:

(о Ч I ч

1 7 7 У 7 "РР~ Xарр , ари = (1 - хг )Б}, аЦ? = (1 - V Щ,

/=пЬ( Р)

app = -

aN =

Sf • dPN

app = -

a

pw _

N

= (1 - Xf)S}, aPpu = 2f

apv =

f=nb( P)

2 If/

apw =

2f

f=nb( P)

f=nb( P)

T

k

k

T

k

k

<

При записи данных коэффициентов использовался алгоритм неортогональной коррекции [10] для возможности правильного расчета на произвольных неструктурированных сетках, а также формула для вычисления давления на грани pf с помощью линейной интерполяции по значениям в центрах ячеек [6]:

Pf = Xfpp + (1 - Х/)рм. Для первого уравнения системы (4) правая часть имеет вид:

bp = Z D f VPf 's f - D f VPf

f=nb( P)

Sf -

S f ' S f

Sf • d

d

PN

PN

(6)

(7)

Суммирование по индексу «/» для второго уравнения системы (4) - уравнения сохранения импульса, дает коэффициенты общей матрицы системы для вычисления компонент скорости:

Г \

ии _ W _ _

UN — им — им —

*N

*N

Za ff

k

S f •S f s f •d

+ min

PN

0, Z a(k )p(k) S k

f

(8)

Первое слагаемое выражений (8) относится к диффузионному слагаемому, и коэффициенты общей матрицы системы будут иметь вид:

aUNP — (1 - Xf )Sf, aN — (1 - Xf )Sy, atf — (1 - V )S}

(9)

auPp —

Z X fSf

f—nb( P)

aj —

Z hS}

f—nb( P)

awp —

Z X fSf

f—nb( P)

Как и для коэффициентов (5) для записи использовался алгоритм неортогональной коррекции.

Второе слагаемое выражений (7) - конвективная составляющая, которая аппроксимируются с помощью любой известной дифференциальной схемы, применимой на произвольных неструктурированных сетках [6, 11, 12]. Во втором выражении уравнений (10) используется противопотоковая схема (Upwind Differences, UD). Также используются: противопоточ-ная схема с линейной интерполяцией (Linear Upwind Differences, LUD), схема QUICK, центрально-разностная схема (Central Differences, CD), схемы семейства NVD (Normalized Variable Diagram), гибридные схемы - все перечисленные схемы, смешанные с противопоточной схемой для увеличения монотонности.

Дискретизация нестационарного слагаемого может быть осуществлена по одной из известных неявных схем [13] - схема Эйлера, схема Адамса-Бэшфорта. Вклад диффузионного и конвективного слагаемых уравнений сохранения импульса относится и к диагональным коэффициентам общей матрицы системы, которые с учетом нестационарного слагаемого дискретизованного с помощью схемы Эйлера будут иметь вид:

ии Up

—- Z

N—NB( P)

ии aN

a

P

+ Z aPk )p P) -, k T

(k )n(k)—

— - Z aN +Z ap )pPk)-

N—NB(P) k

aP

Z„ww aN

N—NB(P)

+ Z aPk )pPk) -. k T

Для второго уравнения системы (4) правая часть имеет вид:

vv

bP = I

f=nb( P)

bP = I

f=nb( P)

bP = I

f =nb( P)

Vu

• ^ f •

Vv

• ^ f •

Vw'

^ f

S, -

S, -

S f -

S f •S f

S f • d PN

S f •S f d

S f • d PN

S f •S f

S f • d PN

d

PN

+

I а ? ^ )Mp- +I а(к )p(k) gxV.

V

(k)n(k),

d

PN

+

I aPk)pPk)VPV +I а(k )p(k) gyV, (П)

d

PN

+

I а (k )p (k) WpV +I а(k )p(k) gzV.

V

(k) (k)

aPP aPu UP apv aPw Up ' Pp UPP aN a Pu aN UPv aN aPW aN PN bP

auPp uu Up apv apw Up + I N=NB( P) UUP aN uu aN uv aN uw aN UN h u bP

af af af afw vP UvP aN vu aN nvv aN vw aN vN bP

aw apu afV „ww Up Wp UWP aN wu aN wv aN ww aN _ _wN _ h w bp

Таким образом, совмещенная система линейных алгебраических уравнений полностью неявного алгоритма для моделирования многофазного потока имеет вид:

(12)

Данная система записана для вычисления общей скорости и общего давления многофазного потока, однако может быть обобщена и на случай расчета, где каждая фаза обладает своей скоростью и физическими свойствами, такими как сжимаемость и турбулентность. Эти обобщения будут проведены в следующих работах.

Для моделирования границ раздела фаз, после решения системы (12), решается уравнение переноса объемных долей (третье уравнение системы (1)), которое решается для (n - 1) объемной доли фаз. Его дискретизация методом конечных объемов осуществляется по схеме полностью аналогичной той, что используется для уравнения сохранения импульса. Для аппроксимации конвективного слагаемого уравнения переноса объемных долей используется схема M-CICSAM [3], относящаяся к классу сжимающих схем высокого разрешения и обеспечивающая сохранение минимально возможной толщины границы раздела сред, а также сохранение формы распределения объемных долей при параллельном переносе и вращении. Дискретизация нестационарного слагаемого также осуществляется по схемам Эйлера или Адамса-Бэшфорта.

В алгебраической форме данная система уравнений для k -й фазы имеет вид:

af ) + I aN 4k} = bP}. (13)

N=NB( P)

Коэффициенты матрицы неявного решения уравнения (13) имеют вид:

aNk) = min (о, u(fl)Sf ), af) =- I a(k} + V,

N=NB( P) T

(14)

b(k) = u(n-l) S (а(к) а(к) )+ V а(к)n-l bP =-uf Sf -(амс - аГО)+_aP

где а^, а^ - значение объемной доли на грани найденной по схеме M-CICSAM и по про-

(к ),н-1

тивопоточной схеме соответственно, ар - значение объемной доли на прошлый шаг по

времени. Представленные слагаемые получены с помощью дискретизации методом отложенной коррекции [10].

Для численного решения итоговая система уравнений должна быть дополнена начальными и граничными условиями. На твердых стенках градиент давления и объемных

др да^

долей равен нулю: — = 0, —к = 0, значение скорости равно нулю: и = 0, V = 0, = 0. На

дп дп

входе задано фиксированное значение скорости и объемных долей, градиент давления равен

k

k

k

k

k

k

нулю: — = 0. На выходе фиксируется статическое давление, градиенты скорости и объем-

дн

ди ду д^ да к

ных долей равны нулю — = 0, — = 0, — = 0, -= 0.

дн дн дн дн

Представленный неявный алгоритм расчета многофазных течений со свободной поверхностью решается с использованием многосеточного метода, позволяющего существенно ускорить вычислительную процедуру [7, 14]. Современный обзор методов ускорения гидродинамических расчетов представлен в [14]. Модель реализована в пакете программ ЛОГОС -российском программном продукте инженерного анализа, предназначенного для решения сопряженных трехмерных задач конвективного тепломассопереноса, аэродинамики и гидродинамики на параллельных ЭВМ [1, 7, 8, 14].

2. Базис задач валидации для течений со свободной поверхностью

Практические все течения со свободной поверхностью нестационарные, поэтому основная цель валидации сводится к необходимости сравнения с экспериментальными данными физических характеристик на определенный момент времени или оценке временного периода их колебаний. Для оценки физических характеристик на определенный момент времени и временного периода их колебаний предлагается использовать следующие задачи валидации, которые мы и назовем «минимальным» базисом.

Таблица 1

Минимальный базис задач валидации

№ Название задачи Оцениваемые параметры Доступные данные

1 Обрушение столба жидкости - Высота столба жидкости по левой стенке резервуара; - положение передней кромки жидкости вдоль дна резервуара - Экспериментальные результаты [20]; - результаты численного моделирования [15]

2 Колебания воды под действием силы тяжести Уровень жидкости на левой стенке резервуара - аналитическое решение [15]; - результаты численного моделирования [15]

3 Обрушение столба жидкости на дно резервуара с препятствием Форма свободной поверхности - Экспериментальные результаты [20]; - результаты численного моделирования [15]

4 Гидравлический удар Высота отраженной волны Результаты численного моделирования [15, 21]

5 Течение через шлюзовые ворота Форма свободной поверхности Результаты численного моделирования [15]

6 Падение шара в жидкость Вертикальная координата центра масс шара - Экспериментальные результаты [17]; - аналитическое решение [17]

7 Падение параллелепипеда в жидкость Амплитуда колебаний уровня воды в бассейне Экспериментальные результаты [17]

8 Падение капли в жидкость Форма свободной поверхности - Экспериментальные результаты [19]; - результаты численного моделирования [18]

3. Описание базиса задач валидации

В данном разделе приводится описание базиса задач валидации, представленных в разделе 2 (табл. 1). Также здесь приводятся результаты валидации, представленной модели, реализованной в пакете программ ЛОГОС.

В расчетах использовалась схема второго порядка точности по времени Адамса-Бэшфорта. Для аппроксимации конвективных слагаемых в уравнении сохранения импульса использовалась противопоточная схема, а для конвективных слагаемых в уравнении переноса объемных долей - схема М-СГСБАМ.

Обрушение столба жидкости

В эксперименте, столб жидкости под действием силы тяжести (воды) падает на дно резервуара [20]. Параметры экспериментальной установки и положение границы раздела фаз (жидкость-воздух) в начальный момент времени приведены на рис. 1 .

открытая граница

а=0Д46 м

воздух

вода

I- 4а -1

Рис. 1. Схема экспериментальной установки

Для проведения численного моделирования была построена структурированная сеточная модель, состоящая из 8400 ячеек. На левую, правую и нижнюю стенки накладывалось граничное условие прилипания, на переднюю и заднюю стенку - скольжение, сверху фиксировалось статическое давление 0 Па.

На рис. 2, а приведены результаты численного моделирования. Сплошной черной линией изображена граница раздела фаз - вода и воздух, черными штрихами показаны вектора скорости. Результаты построены для различных моментов времени: 0,0; 0,2; 0,4; 0,6; 0,8 и 1,0 с. На рис. 2, б показаны результаты численного моделирования, полученные в работе [15]. На рис. 2, в приведены фотографии [15].

В момент времени I = 0,2 с передняя кромка водяного столба встречается с препятствием, формируя мощный всплеск вверх перед препятствием. Высота всплеска, форма образованной каверны хорошо согласуются, как между расчетами, так и экспериментом [15]. Затем образованная всплеском волна достигает правой границы модельной (экспериментальной) области, гидравлический удар об эту стену порождает волну. На момент времени I = 0,6 с приходится самая активная стадия обрушения, образованной волны. При этом расчетные результаты показывают образование воздушного кармана около правой границы области, по форме примерно совпадающего для обоих расчетных случаев. Фотография эксперимента позволяет его идентифицировать. В момент 1,0 с жидкость встречается с левой стенкой.

Из представленных результатов видно, что численные результаты, как по нестационарным характеристикам (время образования и обрушения волны), так и по форме свободной поверхности дают хорошее качественное согласование с экспериментом.

На рис. 3, а представлено изменение высоты столба жидкости вдоль левой стенки резервуара с течением времени, на рис. 3, б - изменение положения передней кромки жидкости

вдоль дна резервуара с течением времени. Высота столба жидкости по левой стенке резервуара, положение передней кромки жидкости вдоль дна резервуара (отсчитываются от левого нижнего угла) и время - величины безразмерные.

t = 0,2 с

t = 1,0 с

а) б) в)

Рис. 2. Обрушение столба жидкости:

а - результаты численного моделирования; б - результаты численного моделирования, полученные в работе [15]; в - экспериментальные данные [20]

Рис. 3. Сравнение результатов численных и экспериментальных исследований:

а - зависимость высоты столба жидкости на левой стенке резервуара от времени; б - зависимость положения передней кромки жидкости от времени

Сплошной черной линией показаны кривые, полученные в результате численного моделирования, штриховой линией - результаты численного моделирования, полученные в [15], маркерами отмечены экспериментальные данные [15].

При вычислении высоты столба жидкости погрешность результатов численного моделирования относительно экспериментальных данных не превышает 1%. При определении положения передней кромки жидкости погрешность результатов численного моделирования относительно экспериментальных данных колеблется от 3,4 до 9,6 %. Видно, что численные результаты дают хорошее качественное и количественное согласование с экспериментом [20] и численным решением, приведенными в работе [15].

Колебания воды под действием силы тяжести

В эксперименте исследуются колебания жидкости, налитой в резервуар, под действием силы тяжести. Параметры численной схемы и положение границы раздела фаз (жидкость-воздух) в начальный момент приведены на рис. 4. В начальный момент времени форма свободной поверхности имеет форму половины периода косинуса. Под действием силы тяжести она начинает колебаться. Теоретический период колебаний Т составляет 0,3739 с [15].

Рис. 4. Схема экспериментальной установки

Для проведения численного моделирования была построена структурированная сеточная модель размером 0,1x0,065x0,01 м, состоящая из 33600 ячеек. На левую, правую и нижнюю стенки накладывалось граничное условие прилипания, на переднюю и заднюю стенку - скольжение, сверху фиксировалось статическое давление 0 Па.

На рис. 5 приведены результаты численного моделирования колебания поверхности жидкости под действием силы тяжести. Сплошной черной линией изображена граница раздела фаз - вода и воздух, серым цветом выделена жидкая фаза, черными штрихами показаны вектора скорости. Результаты построены для различных моментов времени, Т - теоретический период колебаний. Под буквой «а» показаны результаты использования метода, описанного в п. 1, под буквой «б» - полученные в работе [15]. В момент времени I = Т/4 граница раздела фаз принимает вид горизонтальной линии, а затем в I = Т/2 - снова имеет форму полупериода косинусоиды, колеблющейся в противофазе. Процесс повторяется во времени.

На рис. 6 показано колебание уровня жидкости вдоль левой стенки резервуара с течением времени. В табл. 2 приведены погрешности численных результатов по времени относительно теоретических колебаний по закону косинуса. Погрешность по времени вычисляется следующим образом:

• 100%,

где — период колебаний, полученный по результатам численного моделирования, ^ — теоретический период Т [15].

1

1

Видно, что численные результаты дают хорошее качественное и количественное согласование с теорией и численным решением, приведенным в работе [15].

г = 0

г = Т/2

а) б)

Рис. 5. Колебания жидкости под действием силы тяжести:

а - результаты численного моделирования по методу в п. 1; б - результаты численного моделирования, полученные в работе [15]

Рис. 6. Колебания уровня жидкости вдоль левой стенки резервуара под действием силы тяжести

Таблица 2

Погрешность полученных результатов численного моделирования

Относительная погрешность по времени, %

Время Результаты численного модели- Результаты численного моделирова-

рования [15] ния по методу в п. 1

2 • Т 0,00 0,70

4 • Т -0,75 0,76

6 • Т 0,04 0,03

Обрушение столба жидкости на дно резервуара с препятствием

В эксперименте исследуется обрушение столба жидкости (воды) на дно резервуара, где располагается препятствие прямоугольного сечения [20]. Параметры экспериментальной установки и положение границы раздела фаз (жидкость-воздух) в начальный момент времени приведены на рис. 7.

открытая граница

воздух а=0,146 м (1=0,024 м

вода . (1 , 2(1

Стенка

■1 1

4 1

Рис. 7. Схема экспериментальной установки

Для проведения численного моделирования была построена структурированная сеточная модель, состоящая из 92000 ячеек. На левую, правую и нижнюю стенки накладывалось граничное условие прилипания, на переднюю и заднюю стенку - скольжение, сверху фиксировалось статическое давление 0 Па.

На рис. 8, а приведены результаты численного моделирования обрушения плотины. Сплошной черной линией изображена граница раздела фаз - вода и воздух, черными штрихами показаны вектора скорости. Результаты построены для различных моментов времени: 0,1; 0,2; 0,3; 0,4; 0,5 и 0,6 с. На рис. 8, б показаны результаты численного моделирования, приведенные в работе [15]. На рис. 8, в приведены фотографии, сделанные в ходе эксперимента [20].

В момент времени I = 0,1 с передняя кромка водяного столба, построенная по результатам численного моделирования, еще не встретилась с препятствием, однако на экспериментальной фотографии виден тонкий слой жидкости, который уже достиг нижней части препятствия. В момент времени I = 0,2 с гидравлический удар о препятствие порождает волну. На момент времени I = 0,5 с приходится самая активная стадия обрушения образованной волны. При этом расчетные результаты показывают образование воздушного кармана около правой границы области по форме примерно совпадающего для обоих расчетных случаев. Фотография эксперимента позволяет его идентифицировать. Видно, что численные результаты дают хорошее качественное согласование с экспериментом.

t = 0,1 с

а) б) в)

Рис. 8. Обрушение столба жидкости на дно резервуара с препятствием:

а - результаты численного моделирования; б - результаты численного моделирования, полученные в

работе [15]; в - экспериментальные данные [20]

Гидравлический удар

Численно исследуется гидравлический удар, который возникает в процессе движения горизонтального слоя жидкости в ограниченной стенками области. Жидкость, двигаясь с начальной скоростью 1 м/с, ударяется о противоположную стенку, происходит ее обрушение и возникновение обратной волны. Параметры численной схемы приведены на рис. 9 [15]. Длина и высота емкости составляют 12.5 и 4.0 м. соответственно.

открытая граница

Рис. 9. Параметры численной схемы

Для моделирования этой задачи была построена структурированная сеточная модель,

состоящая из 6000 ячеек. На левую, правую и нижнюю стенки накладывалось граничное условие прилипания потока, сверху фиксировалось статическое давление 0 Па.

В табл. 3 приведены значения высоты И0 (рис. 9), которую достигает вода в результате обрушения с учетом g = 1 м/с2 и g = 0,4548 м/с2. Высота И0 измерялась вдоль правой стенки в тот момент времени, когда волна, отразившись от правой границы расчетной области, достигала левой, и поверхность жидкости за фронтом волны была примерно параллельна дну емкости. Была рассчитана относительная погрешность численных результатов относительно аналитического значения. Видно хорошее количественное согласование численных и аналитических результатов.

Таблица 3

Сравнение численных результатов

Расчетный случай g, м/с2 к0, м Относительная погрешность, %

Результаты численного моделирования [21] Результаты численного моделирования по методу в п. 1

1 1,0 2,19 2,20 0,46

2 0,4548 2,90 2,93 1,03

t = 3,50 с

t = 4,90 с

1 = 7,94 с

t = 9,88 с

а) б)

Рис. 10. Гидравлический удар во втором расчетном случае для моментов времени 3,50; 4,90; 7,94 и 9,88 с:

а - результаты численного моделирования; б - опубликованные в работе [15]

На рис. 10 приведены результаты численного моделирования обрушения жидкости с помощью: а) представленной в настоящей работе модели; б) опубликованные в [15]. Сплош-

ной черной линией изображена граница раздела фаз - вода и воздух, серым цветом выделена жидкая фаза, черными штрихами показаны вектора скорости. В момент I = 3,50 с отраженная от правой стенки ударная волна начинает двигаться в обратном направлении. В момент I = 4,90 с фронт волны изгибается и начинает формироваться каверна. Затем в моменты времени I = 7,94 с и I = 9,88 с отраженная волна, сформировав несколько каверн, достигает левой стенки.

Течение через шлюзовые ворота

Численно исследуется задача о течении под тонкой перегородкой, разделяющей бассейны с разными уровнями воды. Параметры численной схемы приведены на рис. 11 [15]. В начальный момент времени вода, заполняющая левый бассейн, под действием силы тяжести начинает перетекать в правый бассейн под тонкой перегородкой, имитирующей шлюзовые ворота.

открытая граница

А 1 -ч- ® о L воздух

1 вода

С Стенки

- 0.3 м - - 1.61 Л1 -==-

Рис. 11. Параметры численной схемы

Для моделирования этой задачи была построена структурированная сеточная модель, состоящая из 6000 ячеек. На левую, нижнюю и нижнюю часть правой стенки накладывалось граничное условие прилипания потока, сверху и на верхней части правой стенки фиксировалось статическое давление 0 Па.

На рис. 12 приведены результаты моделирования течения: а) результаты численного моделирования; б) результаты, опубликованные в работе [15]. Сплошной черной линией изображена граница раздела фаз - вода и воздух, черными штрихами показаны вектора скорости.

t = 0,25 с

t = 0,70 с

t = 0,93 с

а) б)

Рис. 12. Течения через шлюзовые ворота:

а - результаты численного моделирования; б - опубликованные в работе [15]

В момент времени I = 0,25 с наблюдается формирование первой каверны сразу за тонкой стенкой. В момент времени I = 0,70 с формируется третья каверна, в I = 0,93 с фронт волны достигает правой стенки.

Следует отметить, что в работе [15] не указано, в какие моменты времени были сделаны фотографии. Авторами этой статьи были самостоятельно выбраны моменты времени для построения результатов численного моделирования. Видно хорошее качественное согласование численных и аналитических результатов.

Падение шара в жидкость

В задаче моделируется падение твердого шара в воду с небольшой начальной скоростью. Подобный эксперимент проводился для шаров из разных материалов и изложен в работе [17].

Параметры численной схемы приведены на рис. 13. Шар радиусом г = 1,27 см непосредственно перед столкновением с водной поверхностью имеет скорость V = 2,17 м/с. Уровень воды к в начальный момент времени составляет 0,2 м. Расчетная область представляет собой цилиндр радиусом Я = 0,25 м и высотой Н = 0,2254 м.

Рассматривалось два расчетных случая: в первом случае шар был сделан из полипропилена и плотность его составляла 0,86 от плотности воды, во втором - из стали с плотностью 7,86 от плотности воды.

Рис. 13. Параметры численной схемы

Для моделирования этой задачи была построена структурированная сеточная модель, измельченная в области движения твердого тела и состоящая из 800000 ячеек (см. рис. 14). На левую, правую и нижнюю стенки накладывалось граничное условие прилипания потока, сверху фиксировалось статическое давление 0 Па.

а)

б)

а - общий вид; б -

Рис. 14. Дискретная модель:

измельчение сетки в области движения твердого тела

На рис. 15 показана динамика погружения стального шара в воду: а) экспериментальные фотографии [17]; б) результаты численного моделирования, построенные для соответствующих моментов времени. Серым цветом выделена жидкая фаза, сплошной черной линией отмечена граница раздела фаз, шар выделен черным цветом. В момент I = 5,9 мс шарик полностью погружается под воду. В момент I = 54,9 мс около над поверхностью шара начинает формироваться каверна, формирование которой к моменту времени I = 68,9 мс заканчивается по результатам числененого моделирования и полностью завершается в эксперименте.

Видно хорошее качественное согласование численных и экспериментальных результатов.

5.9 ms 12.9 19.9 26.9 33.9 40.9 47.9 54.9 61.9 68.9 75.9

б)

Рис. 15. Погружение стального шара в воду:

а - экспериментальные фотографии [17]; б - результаты численного моделирования

На рис. 16 показано изменение глубины погружения шара в воду с течением времени. Приведены результаты для двух расчетных случаев. Теоретические результаты приведены в работе [17].

Рис. 16. Изменение вертикальной координаты центра масс шарика от времени

Погрешность результатов моделирования относительно теоретических результатов составляет:

• для полипропиленового шара - 5,8%;

• для стального шара - 6,2%.

Падение параллелепипеда в жидкость

Верификацию представленного алгоритма на правильное описание распространения волн на свободной поверхности можно осуществить с помощью численного моделирования эксперимента, описанного в [16]. В этом эксперименте с высоты Н вдоль торцевой стенки бассейна в воду свободно падает прямоугольный параллелепипед с высотой Н\ и длиной I. Прямоугольный бассейн с горизонтальным дном заполнен водой на высоту к < Н\. В натурном эксперименте длина бассейна составляла 4,3 м, а ширина 0,2 м. В невозмущенном состоянии вода находится в состоянии покоя. Расчетная модель представлена на рис. 17.

1

н ? ? 1 1 1 Л 1 1

Д-„ = 16 ЛГ„ = 31,5 *0,м Рис. 17. Параметры численной схемы

Рассматривалось несколько вариантов постановок задачи, отличающихся друг от друга размерами тела, высотой, с которой оно падает и глубиной заполнения бассейна. В безразмерной форме эти параметры выглядят следующим образом [16]:

№ = H/h, Hi0 = Hi/h, f = l/h, p0 = (pi - p) / p, X = x/h,

где p1, p - плотность твердого тела и плотность жидкости соответственно. Величины h, H0,

0 0 0 0

Hi l , p , x меняются в зависимости от постановки.

В первом случае значения этих величин следующие: h = 8 см, H° = 3,75, Hi0 = 2,26, l0 = 1,15, p0 = 0,215.

Для второго расчетного случая h = 8 см, H0 = 2,90, Hi° = 2,26, l0 = 0,575, p0 = 0,215. На расстоянии x0 = 16 установлен мареограф - датчик, фиксирующий колебания водной поверхности.

Для третьего расчетного случая h = 4 см, H° = 4,75, H10 = 4,5, f = 1,15, p0 = 0,215, мареограф был установлен на расстоянии x0 = 31,5.

Для моделирования падения параллелепипеда в воду была построена структурированная расчетная сетка, состоящая примерно из 1 млн ячеек в области с размерами 4,3x0,4x0,5 м. В месте падения параллелепипеда и установки мареографа сетка более мелкая. На левую, правую и нижнюю стенки дискретной модели накладывалось граничное условие прилипания потока, сверху фиксировалось статическое давление 0 Па.

На рис. 18 показано образование волн от падения параллелепипеда для первой постановки задачи. Черной линией обозначен контур параллелепипеда. Серым цветом выделена жидкая фаза. В момент времени t = 0,25 с наблюдается ранняя стадия развития каверны и брызговой струи. В момент времени t = 0,3 с струя продолжает развиваться. Через 0,1 с дно канала обнажилось, скорость формирования струи уменьшилась. В момент t = 7,5 с показана заключительная стадия схлопывания каверны и генерации брызговой струи. Следует

отметить, что в работе [16] не указано, в какие моменты времени были сделаны фотографии. Авторами этой статьи были самостоятельно выбраны моменты времени для построения результатов численного моделирования.

б)

Рис. 18. Образование волн от падения параллелепипеда для первой постановки задачи:

а - результаты численного моделирования, построенные в моменты времени 0,1; 0,2; 0,3 и 0,4 с;

б - экспериментальные данные [16]

0,06

-0,02

Время, с

Рис. 19. Амплитуда колебаний уровня воды в бассейне для второй постановки задачи

г <

-0,02

Время, с

Рис. 20. Амплитуда колебаний уровня воды в бассейне для третьей постановки задачи

На рис. 19 и рис. 20 приведены зависимости амплитуды колебаний уровня воды в бассейне от времени для второй и третьей постановки задачи соответственно. Точками на гра-

фик нанесены экспериментальные данные, записанных неподвижным волномером [16]. Сплошной черной линией показана амплитуда колебаний уровня воды, налитой в бассейн, полученных по результатам численного моделирования.

Сравнение мареографных записей показывает хорошее совпадение для первой приходящей волны. Вычисленная амплитуда вторая волны («отрицательной») также хорошо соответствует экспериментальным данным, однако наблюдается различие для дальнейших колебаний поверхности, которые имеют миллиметровую амплитуду. Это может быть связано как с качеством расчетной сетки, так и со способом записи численных данных в мареографной точке.

Падение капли в жидкость

Рассматривается задача об ударе капли о поверхность жидкости с образованием короны. Капля имеет те же параметры, что и жидкость: плотность 1200 кг/м , вязкость 0,022 Н-с/м2, коэффициент поверхностного натяжения 0,0652 Н/м, числа Вебера и Рейнольд-са в момент столкновения капли с поверхностью воды равны 2010 и 1168 соответственно. Радиус капли составляет 0,0021 м. Результаты экспериментальных и численных исследований приведены в работах [18] и [19] соответственно. Расчетная область представляет собой цилиндр с высотой Н 0,15 м и радиусом Я 0,03 м (см. рис. 21).

открытая граница

Рис. 21. Численная схема

1 = 0,25 с

t = 1,0 с

t = 3,0 с

а) б) в)

Рис. 22. Динамика удара капли о поверхность жидкости

Для моделирования этой задачи была построена структурированная сеточная модель с измельчением в области падения капли и образования короны, состоящая из 2,8 млн ячеек. Сверху фиксировалось статическое давление 0 Па.

На рис. 22 приведена динамика удара капли о поверхность жидкости: а) и б) результаты численного моделирования падения капли; в) экспериментальные фотографии [19]. Сплошной черной линией изображена граница раздела фаз - вода и воздух, серым цветом выделена жидкая фаза. В момент t = 0,25 с капля начинает погражуться в воду. В момент t = 1,0 с видна корона брызг, сформировавшаяс в результате падения капли. В момент времени t = 3,0 с корона продалжает увеличиваться.

Видно хорошее качественное согласование численных и экспериментальных результатов.

Заключение

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

С использование сформированного базиса задач была проведена валидация пакета программ ЛОГОС для данного класса задач. Относительная погрешность результатов не превышает 10%. Видно хорошее качественное и количественное согласование численных и экспериментальных результатов. По результатам исследований можно сделать вывод, что ЛОГОС позволяет с достаточной точностью моделировать течения со свободной поверхностью.

Библиографический список

1. Козелков, А.С. Минимальный базис задач для валидации методов численного моделирования турбулентных течений вязкой несжимаемой жидкости / А.С. Козелков [и др.]/ / Труды Нижегородского государственного технического университета им. Р.Е. Алексеева. 2014. № 4(104). С. 21-69.

2. AIAA, «Guide for the Verification and Validation of Computational Fluid Dynamics Simulations», AIAA G-077-1998.

3. Waclawczyk, T. Remarks on prediction of wave drag using VOF method with interface capturing аpproach / T. Waclawczyk, T. Koronowicz // Archives of civil and mechanical engineering. 2008. V. 8. P. 5-14.

4. Волков, К.Н. Течения газа с частицами / К.Н.Волков, В.Н. Емельянов. - М.: Физматлит, 2008. - 600 с.

5. Kolev, N.I. Multiphase Flow Dynamics / N.I. Kolev. - Berlin, Heidelberg: Springer-Verlag, 2007.

6. Ferziger, J.H. Computational methods fluid dynamics / J.H. Ferziger, M. Peric. - Berlin, Heidelberg: Springer, 2001.

7. Козелков, А.С. Реализация метода расчета вязкой несжимаемой жидкости с использованием многосеточного метода на основе алгоритма SIMPLE в пакете программ ЛОГОС / А.С. Козелков [и др.]/ // ВАНТ, сер. Математическое моделирование физических процессов. 2013. Вып. 4.

8. Darwish, M. A coupled finite volume solver for the solution of incompressible flows on unstructured grids / M. Darwish, I. Sraj, F. Moukalled // Journal of Computational Physics. 2009. V. 228. P.180-201.

9. Rhie, C.M. A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation / C.M. Rhie, / W.L. Chow // AIAA. 1983. V. 21. P. 1525-1532.

10.Jasak, H. Error Analysis and Estimation for the finite volume method with applications to fluid flows. Thesis submitted for the degree of doctor // Department of Mechanical Engineering, Imperial College of Science, 1996.

11. Математические модели и алгоритмы для численного моделирования задач гидродинамики и аэродинамики: учеб. пособие / А.С. Козелков [и др.]. - Нижний Новгород: НГТУ, 2014. - 166 с.

12.Козелков, А.С. Моделирование течений вязкой несжимаемой жидкости на неструктурированных сетках методом отсоединенных вихрей / А.С. Козелков [и др.] // Математическое моделирование. 2013. Т. 26. № 8. C. 81-96.

13.Роуч, П. Вычислительная гидродинамика / П. Роуч. - М.: Мир, 1980. - 616 с.

14.Методы ускорения газодинамических расчетов на неструктурированных сетках / А.С. Козелков [и др.]. - М.: Физматлит, 2013.

15.Ubbink, O. Numerical prediction of two fluid systems with sharp interfaces. Thesis. Imperial College of Science, Technology & Medicine. London. 1997.

16.Букреев, В.И. Гравитационные волны при падении тела на мелкую воду / В.И. Букреев, А.В. Гусев // Прикладная механика и техническая физика. 1996. Т. 37. №2. С. 90-98.

17.Aristoff, J.M. The water entry of decelerating spheres / J.M. Aristoff [et al.] // Phys. Fluids. Am. Inst. Phys. 2010. No. 22.

18.Джанг, Т. Численное исследование динамики удара капли о поверхность жидкости с образованием короны / Т. Джанг, Д. Оуянг, Х. Ли, Д. Рен, С. Ванг // Прикладная механика и техническая физика. 2013. Т.54. № 5. С. 38-47.

19.Wang, A.B. Splashing impact of a single drop onto very film liquid films / A.B. Wang, C.C. Chen // Phisics of Fluids. 2000. V. 12. No. 9. P. 2155-2158.

20.Koshizuka, S. A particle method for incompressible viscous flow with fluid fragmentation / S. Ko-shizuka, H. Tomako, Y. Oka // Computational Fluid Dynamics Journal. 1995. V. 4(1). P. 29-46.

21.Hirt, C.W. Volume of Fluid Method (VOF) for the dynamics of free boundaries / C.W. Hirt, B.D. Nichols // Journal of Computational Physics. 1981 V. 39. P. 201-225.

Дата поступления в редакцию 06.05.2015

A.S. Kozelkov1'2, A.A. Kurkin2, I.L. Sharipova1, V.V. Kurulin1, E.N. Pelinovsky2, E.S. Tyatyushkina1, D.P. Meleshkina1, S.V. Lashkin1, N.V. Tarasova1

MINIMAL BASIS TASKS FOR VALIDATION OF METHODS OF CALCULATION

OF FLOWS WITH FREE SURFACES

FSUE «RFNC - VNIIEF»1, Nizhny Novgorod state technical university n.a. R.E. A^eev

Purpose: In this paper the validation process, which is an important step towards the industrial application of engineering software packages is discussed.

Method: The study is based on physico-mathematical and numerical models designed to simulate the flow of a viscous incompressible fluid with free surfaces.

Results: In order to improve and systematize the knowledge in this paper a basis of validation tasks required to assess the accuracy of the simulation software solutions designed for the simulation of flows of incompressible viscous fluid with free surfaces are developed. Held its systematization and generalization.

Application domain: Presented results allow the developer of engineering software packages focus on the computation of error modeling, rather than searching for reliable data.

Key words: validation, engineering software packages, flow, free surface, numerical simulation.