CC BY
77
Вычислительные технологии
Том 18, № 3, 2013
Математическое моделирование задачи прорыва дамбы
А. Ж. ЖАйНАКОВ1, А. Ы. КУРБАНАЛИЕВ2 1 Институт горного дела и горных технологий им. У. И. Асаналиева Кыргызского государственного технического университета им. И. Раззакова,
Бишкек, Кыргызстан 2Кызыл-Кийский гуманитарно-педагогический институт, Баткенского государственного университета, Кызыл-Кия, Кыргызстан e-mail: jainakov-41@mail.ru, kurbanaliev@rambler.ru
Представлены результаты математического моделирования классической задачи прорыва дамбы. Основу трёхмерного нестационарного моделирования составляют осреднённые по Рейнольдсу уравнения Навье — Стокса. Для отслеживания положения свободной границы применяется метод объёма жидкости. Адекватность математической модели проверена путём сравнения с экспериментальными данными.
Ключевые слова: прорыв дамбы, свободная граница, метод объёма жидкости, ОрепРОЛМ.
Введение
Гидротехнические сооружения, относящиеся к числу сложных технологических объектов, создают целый ряд социальных, экологических и экономических проблем даже при нормальном режиме работы. В результате разрушения таких сооружений как плотины или дамбы происходят значительные затопления местностей, расположенных ниже по течению реки, с находящимися на них населёнными пунктами и разного рода хозяйственными объектами.
В настоящее время для определения параметров прорывных волн применяются различные методы, например, натурные исследования, физический эксперимент, численное моделирование в одномерной или двумерной постановке задачи. В пределах своих границ применимости они позволяют предсказать время появления переднего фронта волны, границы зон затопления, глубину и продолжительность затопления местности. При этом актуальной задачей является определение гидродинамического воздействия переднего фронта волны прорыва на различные инженерные сооружения. Вместе с тем существующие методы численных расчётов и модели тех или иных физических явлений нуждаются в совершенствовании.
Предлагаемая в настоящей работе математическая модель основана на усреднённых по Рейнольдсу трёхмерных нестационарных уравнениях Навье — Стокса, замыкание которых осуществляется при помощи стандартной к — е-модели турбулентности. Для определения положения изменяющейся с течением времени свободной
поверхности применялся метод объёма жидкости, имплементированный в многофазный решатель т1егЕоат открытого пакета ОрепЕОАМ 1.7.1. [1].
1. Математическая модель
Законы сохранения массы и импульса для несжимаемой вязкой изотермической жидкости при отсутствии массовых и поверхностных сил приводят к следующим, усреднённым по Рейнольдсу нестационарным уравнениям Навье — Стокса [2, с. 293]:
д
дХ (рЩ) = 0, (1)
д / . д ( _ . "^Л дР . дТг
ж(рщ) + дх +рщи^ = - дх +
1]
где и; — ¿-я компонента средней скорости в направлении координаты хг, р — плот-
_ _ дщг дЩ
ность, р — среднее давление, т; = и I —--+ —— I — средний тензор вязких напря-
1 дх; дхг
жений, и — молекулярная динамическая вязкость, рп'ги] — напряжения Рейнольдса, требующие моделирования. Осреднение производится по времени, а штрих означает флуктуационную часть. При наличии массовых и других сил необходимо уравнение (2) дополнить соответствующими членами.
Для замыкания систем уравнений (1)-(2) необходимо ввести ту или иную модель турбулентности. Многие из моделей турбулентности, используемые на практике, основаны на понятиях турбулентной вязкости и турбулентной диффузии. Для течений общего вида введенная Буссинеском турбулентная вязкость, связывающая напряжения Рейнольдса с градиентами осреднённого течения, может быть записана в форме [2, с. 294]
дЩг дщ\ 2
-рЩЩ = М дх+дх - 3рк6г;,
х ч
к2
где и = рСп--турбулентная вязкость.
£
Кинетическая энергия турбулентности к и скорость ее диссипации £ определяются из транспортных уравнений [2, с. 295-296]
д(рк) д(рЩк) д ( иЛ дк с
+--о-= ТГ- и +--й--+ Рк - р£, (3)
дЪ дх ; дх ; \ ак) дх
д(р£) , д(рЩ£) _ ^ £р пг £\ д (иЛ д£
+--^- = Се!Т Рк — РСе2~Г + Т;- -
л , 1 г\ ^ Ы 7 к г Ь2 7 1 г\ I \/
дъ дх ; к к дх ; \iJeJ дх ;
(&Щ дщ \ дЩг
—--+ —— ---скорость генерации энергии турбулентности средним
дх дхг дх
течением. Коэффициенты модели имеют следующие стандартные значения: См = 0.09, Се1 = 1.44, Се2 = 1.92, ак = 1.0, а£ = 1.3.
1.1. Метод объёма жидкости
При моделировании рассматриваемого класса течений особое место занимает метод определения границы раздела между двумя фазами — водой и воздухом. Согласно основной идее метода объёма жидкости [3], для каждой вычислительной ячейки находится некоторая скалярная величина, представляющая собой степень заполнения данной ячейки одной фазой, например водой. Если в какой-то ячейке эта величина равна 0, то она пустая, если равна 1, то заполнена, если её значение лежит между 0 и 1, то можно сказать, что ячейка содержит свободную (межфазную) границу. Другими словами, объёмная доля воды а определяется как отношение объёма воды в ячейке к полному объёму данной ячейки. Соответственно величина 1 — а представляет собой объёмную долю второй фазы в данной ячейке — воздуха. В начальный момент времени Ь = 0 даётся распределение поля этой величины а{х,у,г,Ь = 0) (см. ниже рис. 2), и дальнейшая её временная эволюция вычисляется как решение транспортного уравнения
да д {аЩ)
Ж + ~дх- = 0 (5)
Положение свободной границы определяется уравнением а{х,у, г,Ь) = 0. Поэтому физические свойства газожидкостной смеси находятся усреднением с соответствующим весовым коэффициентом [2, с. 383]:
р = ар\ + {1 — а)р2, ^ = а/л\ + {1 — а)^2 ■ Здесь индексы 1 и 2 означают жидкую и газовую фазы.
1.2. Начальные условия
Для нестационарной задачи необходимо задание начальных значений всех зависимых переменных. Значения всех компонент скорости равны нулю, так как по условию рассматриваемой задачи до момента времени Ь = 0 движение отсутствует. Давление тоже равно нулю. Кинетическая энергия турбулентности и скорость её диссипации имеют некоторые малые значения, что обеспечивает хорошую сходимость численного решения на первых шагах интегрирования. Начальное распределение объёмной доли а неоднородно, поскольку часть расчётных ячеек была заполнена водой (см. ниже рис. 2).
1.3. Граничные условия
На неподвижных твёрдых стенках расчётной области задано условие прилипания, что определяет равенство нулю всех компонентов вектора скорости. Для давления и объёмной доли воды заданы условия непроницаемости — нулевой градиент этих величин (условие Неймана — Ур = 0, У а = 0); для кинетической энергии турбулентности к и скорости её диссипации е граничные условия задавались с помощью аппарата пристеночных функций [2, с. 298-299].
На границе раздела двух фаз необходимо задание граничных условий для давления и скорости в следующей форме [4, с. 6]:
—р + 2^-
ди„
дп
—ро + аК, ^
дип дщ дЬ дп
0,
где un,ut — нормальная и тангенциальная компоненты скорости соответственно, p0 — атмосферное давление, a — коэффициент поверхностного натяжения, K — кривизна поверхности.
2. Численная модель
Дискретизация расчётной области выполняется методом контрольных объёмов [5, с. 24]. Основная идея этого метода состоит в том, что вся расчётная область разбивается на ячейки гексаэдральной формы с последующей заменой дифференциальных уравнений алгебраическими, представляющими собой интегральные балансовые соотношения для каждой ячейки.
Для вычисления интегралов по контрольному объёму использовался метод Гаусса, а соответствующие значения величин на гранях контрольного объёма вычислялись из значений в центрах соседних ячеек путём применения тех или иных интерполяционных схем. Сведения о методах дискретизации и решения алгебраических уравнений приведены в табл. 1, 2. Более детальная информация об обсуждаемом вопросе дана в [1].
Для решения уравнения переноса (5) для объёмной доли жидкости а применяется созданный в компании OpenCFD решатель MULES, использующий многомерный универсальный ограничитель для явного решения и обеспечивающий ограниченность объёмной доли фазы, независимой от применяемой основной численной схемы, структуры сетки и т. д.
Таблица 1. Схемы дискретизации
Уравнение Производная по Конвективный Диффузионный
времени член член
Момент импульса Euler Gauss linear Gauss linear
Объёмная доля жидкости, а То же Gauss vanLeer То же
Кинетическая энергия тур-
булентности, к — » — Gauss upwind — » —
Скорость диссипации энер-
гии, £ — » — То же — » —
Таблица 2. Решатели уравнений для основных переменных
Уравнение Решатель Точность
Коррекция давления Метод сопряжённых градиентов
с предобуславливанием PCG 1e-10
Давление То же 1e-7
Скорость Метод бисопряжённых градиентов
с предобуславливанием PBiCG 1e-6
Кинетическая энергия тур-
булентности, к То же То же
Скорость диссипации энер-
гии, £ — » — — » —
С целью проверки независимости численных расчётов от размера сетки расчёты проводились на разных сетках: 60 х 25 х 25, 90 х 40 х 40, 135 х 60 х 60. Представленные в работе результаты соответствуют последней сетке. Отдельные расчёты были проведены также для сеток 210 х 110 х 20 и 300 х 210 х 20.
3. Результаты численных расчётов
Геометрия модельной задачи и принятая система координат представлены на рис. 1. Вертикальная координатная ось г перпендикулярна к плоскости хОу. Модель представляет собой открытый сверху резервуар в форме гексаэдра длиной 3.22 м, высотой 1.2 м, шириной 1 м [6]. В левом нижнем углу расположен столб воды высотой Н = 0.6 м, длиной 1.2 м, шириной 1 м.
Для измерения силы давления потока жидкости на правую стенку резервуара соответствующие измерительные датчики были расположены в точках Р1, Р2, Р3, Р4, Р5, Р6 и Р7, координаты которых приведены в табл. 3.
Кроме того, в четырёх сечениях при XI = 2.725, х2 = 2.228, х3 = 1.73 и х4 = 0.6 м были измерены уровни свободной поверхности воды Н1, Н2, Н3 и Н4 соответственно.
На рис. 2 представлены результаты численных расчётов объёмной доли воды для различных моментов времени. При £ = 0 перегородка мгновенно убирается и столб воды под действием силы тяжести устремляется в правую пустую часть резервуара. В момент времени около £ = 0.65 с вода достигает правой стенки и, ударяясь о неё под действием силы инерции, устремляется вверх.
По мере движения вверх по правой стенке поток утончается, и в момент времени £ = 1.4 с, когда сила тяжести превышает силу инерции, происходит обратное движение воды с образованием характерного крючка. Сформированная таким образом обратная волна, достигая основной поток и ударяясь о него, образует вторичную
НА 1.20 НЪ Н2 н\
Рис. 1. Геометрия тестовой задачи
Таблица 3. Координаты измерительных датчиков
Координаты Р1 Р2 Р3 Р4 Р5 Р6 Р7
X 3.22 3.22 3.22 3.22 3.22 3.22 3.22
У 0.160 0.160 0.584 1.000 0.160 0.160 0.160
г 0.164 0.500 0.500 0.500 0.355 0.645 0.790
Рис. 2. Движение столба воды в различные моменты времени
волну и т.д. Примерно около Ъ = 3.0 с обратная волна достигает левой стенки и, отражаясь от неё, опять двигается в сторону правой стенки. После момента времени ъ = 4. 0 с инерция воды значительно уменьшается и дальнейшее рассмотрение движения не представляет практического интереса.
На рис. 3 представлены численные (сплошные линии) и соответствующие экспериментальные (пунктир) данные по высоте водного столба в сечениях с координатами х1 = 2.725м (а), х2 = 2.228м (б), х3 = 1.73м (в) и х4 = 0.6м (г). Видно, что совпадение между приведенными расчётными и экспериментальными данными достаточно удовлетворительно. В целом характер изменения численных данных соответствует эксперименту [6]. Вместе с тем максимальные расчётные значения в точках х1 и х2 для момента времени около ъ = 2 с оказались несколько завышены. В этот момент отражённая от правой стенки волна, сталкиваясь с набегающим течением, приводит к сложному взаимодействию жидкости и воздуха, в результате чего возникает газожидкостная смесь (см. рис. 2, момент времени Ъ = 1.9 с).
Другая причина определяется трудностью воспроизведения экспериментальных условий при численном моделировании. Как указано в [4, с. 19], в эксперименте [6] перегородка удалялась со скоростью 4м/с, а при численном моделировании предпо-
лагается, что она убирается мгновенно. Кроме того, из рис. 3, г следует, что начальная высота водного столба немного меньше 0.6 м. Это может быть причиной того, что численные значения высоты водного столба в точках х1, х2, хз и х4 несколько выше соответствующих экспериментальных величин (см. рис. 3).
Рис. 3. Высота водного столба при различных сечениях
Рис. 4. Давление потока на правую стенку резервуара
Сравнение численных результатов (сплошные линии) и эксперимента (пунктир) по давлению потока на правую стенку резервуара приведено на рис. 4. Видно, что достаточно хорошее совпадение данных имеет место в точке Р3. При численном моделировании имеющийся в эксперименте первый пик не обнаружен, что, вероятно, связано с особенностью переднего фронта потока воды прорыва. Подобные результаты получены также в работах [4, 7].
На рис. 5, взятом из [4], представлены результаты численных расчётов, выполненных с помощью программы СошПо [8], и соответствующие экспериментальные данные для точек Р1, Р2, Р3, Р5, Р6 и Р7. Из рис. 4 и 5 следует, что предлагаемая в настоящей работе математическая модель в сравнении с расчётами [4] лучше описывает рассматриваемый класс течений.
Рис. 5. Сравнение численных и экспериментальных данных [4]
Таким образом, результаты численного моделирования классической задачи прорыва дамбы на основе нестационарных трёхмерных уравнений Навье — Стокса, описывающих динамику газожидкостной смеси со свободной границей, показали хорошее совпадение с данными эксперимента. В дальнейших исследованиях, по мнению авторов, особое внимание необходимо уделить точности расчёта поля давления непосредственно перед ударом потока воды о правую стенку.
Список литературы
[1] HTTp://www.openfoam.org/archive/1.7.1/docs/. OpenFOAM 1.7.1. User Guide.
[2] Ferziger J.H., Peric M. Computational Methods for Fluid Dynamics. Springer-Verlag, 2002. 423 p.
[3] Hirt C., Nichols B. Volume of fluid (VOF) method for the dynamics of free boundaries // J. of Comput. Phys. 1981. Vol. 39. P. 201-225.
[4] Fekken G. Numerical Simulation of Greenwater Loading on the Foredeck of a Ship. MSc. Thesis. Univ. of Groningen, 2004. Available at: http://www.math.rug.nl/~veldman/ preprints/MSc-Fekken.pdf
[5] Патанкар С. Численные методы решения задач теплообмена и динамики жидкости: Пер. с англ. М.: Энергоатомиздат, 1984. 152. c.
[6] Buchner B., van Ballegoyen G. Joint Industry Project: F(P)SO Green Water Loading. MARIN. May 1997. Vol. A3: Scale Effect Tests.
[7] Zhou Z.Q., Kat J.O.D., Buchner B. A nonlinear 3-D approach to simulate green water dynamics on deck // 7th Intern. Conf. on Numer. Ship Hydrodynamics. Nantes, 1999.
[8] HTTp://www.math.rug.nl/~veldman/comflo/comflo.html
Поступила в 'редакцию 6 сентября 2012 г., с доработки — 25 апреля 2013 г.
