CC BY
7УДК 626.421.4
И. В. Липатов, к. т. н., доцент, ВГАВТ.
603600, Нижний Новгород, ул. Нестерова, 5,1_Lipatov@mail.ru.
ФОРМУЛИРОВКА ПОЛНОЙ ПРОСТРАНСТВЕННОЙ МОДЕЛИ ДИНАМИКИ РАЗВИТИЯ ВОЛНОВОГО ПРОЦЕССА В ВЕРХНЕМ ПОДХОДНОМ КАНАЛЕ ШЛЮЗА
В статье формулируется полная трехмерная математическая модель развития волнового процесса в верхнем подходном канале шлюза. Модель учитывает динамику изменения свободной поверхности и турбулентные процессы в потоке. Для модели сформулированы граничные и начальные условия, модель турбулентности. В конце представлены соображения по выбору программного обеспечения и необходимом вычислительном ресурсе для полноценной практической реализации модели.
На сегодняшний день шлюзы стали неотьемлимой частью каждого крупного гидроузла на равнинной реке. С одной стороны они стали способствовать интенсификации развития водный путей, а с другой стали тормозящим фактором при движении судов, находящихся на акватории шлюза и источником повышенной аварийности. Одной из первопричин как первого так и второго недостатков являются сложные, нестационарные волновые процессы развивающиеся в подходных каналах в непосредственной близости от рабочих ворот шлюзов при наполнении-опорожнении камеры.
Естественным путем улучшения ситуации с организацией.судоходства на пришлюзовых участках является изучение гидродинамических процессов как таковых. А затем уже на базе этих знаний строить максимально приемлемую эксплуатационную схему. В данной статье пред по латается остановиться на реализации наиболее современного подхода к исследованию волновых процессов в верхних подходных каналов шлюзов - численному моделированию этих процессов.
Как показал анализ существующих работ в этом направлении, традиционным путем решения подобных задач является интегрирование системы уравнений Сен-Венана. В основе последней лежит предположение о однонаправленности волны в исследуемой области. Такой допуск вполне приемлем, если сечения, для которых необходимо рассчитать изменения глубин находится на значительном расстоянии от первоисточника возмущения. Это класс волновых задач при решении которых фактор влияния граничных условий на конечный результат в интересующей области интегрирования не существенен. То есть влиянием изменения характера эпюры скоростей по глубине потока можно пренебречь. Такую ситуацию в гидротехнике мы имеем при расчете волн паводка, когда моделируем участок в несколько километров. Погрешность в задании граничных условий на входе локализуется в пределах 10-20 размеров глубины потока и существенно не искажает конечную картину волнового процесса в сечениях отстоящих от входного сечения на расстоянии 100 метров и более.
В случае же с волнами формирующимися при поднятии рабочих ворот шлюза, особенно в непосредственной близости к ним, ситуация становится качественно иной. В начальный момент маневрирования отверстие вытекания находится на самом дне канала, так как степень открытия ворот не велика. В середине и конце периода маневрирования, когда сечение для вытекания потока становится соизмеримо с вертикальным размером течения по глубине, и уже идет влияние более чем на 40 % глубины потока, что существенно корректирует ход процесса волнообразования. Амплитуда раскачки волн в начальном сечении становится значительнее, а волнообразования интенсивнее. Таким образом, усредняя скорость на всем сечении потока по глубине,
как это предполагает уравнение Сен-Венана, мы уже на постановочном этапе решения искажаем волновую картину в сторону сглаживания острых углов. Особенно это неприемлемо для участков непосредственно примыкающих к воротам, так как согласно принципу Сен-Венана - ошибка в задании граничных условий постепенно сходится на нет с удалением от места граничной области. С таким явлением мы и имеем дело при решении задачи о волнах паводка.
В случае же исследования гидродинамических процессов в непосредственной близости к воротам (что и вызывает основной эксплуатационный интерес) решение на базе уравнений Сен-Венана становится некорректным по определению. Таким образом возникает необходимость в использовании полной системы гидродинамических уравнений Навье-Стокса:
ГаЧ дг\х г 1 Эр.
V [ахг 3/ Эг2 ^ " р Зх ’
^ V (д\ д*\г г 1 Эр .
а [ах2 Зу2 - ГУ У Р ду'
ауг (дгу1 а2уг д2уг
<31 У\ Эх2 + ду2 + дгг
Зх ду дг р Л ’
Принимая во внимания тот факт, что моделирование процесса волнообразования связано с математическим моделированием поведения двухфазной среды (вода-воздух), для уравнений (1) используем дискретизацию с помощью метода контрольного объёма / 4 /. Суть метода вкратце состоит в следующем. Расчётная область разбивается на некоторое число непересекающихся контрольных объёмов таким образом, что каждая узловая точка содержится в одном контрольном объёме. Дифференциальные уравнения (1) интегрируют по каждому контрольному объёму. Для вычисления интегралов используют кусочные профили, которые описывают изменение исследуемой величины между узловыми точками. В результате находят дискретный аналог дифференциального уравнения, в который входят значения исследуемой величины в нескольких узловых точках. Одним из важных преимуществ метода контрольного объёма является то, что в нём заложено точное интегральное сохранение таких величин как: масса, количество движения и энергия для любого объёма, группы объёмов и всей расчётной области. Это свойство справедливо при любом количестве узловых точек, а не только для предельного случая их очень большого числа. Поэтому даже решение на грубой сетке будет удовлетворять точным интегральным балансам, Принимая во внимание то, что при решении нам необходимо будет отслеживать границы раздела сред вода-воздух, систему уравнений (1) дополним еще одним выражением для пассивного скаляра УОР. С помощью этой формальной величины, в последствии предполагается осуществлять объемное отслеживание распределения сред по расчетной области. Для создания модели предполагается применить метод статической сетки объемного слежения. Первоначальная методика объемного слежения - метод УОР была предложена в работе Херлоу /11/, согласно ей отдельные частицы слежения используются в каждой ячейке. Этот постулат аналогичен принципу метода контрольного объема, где любая переменная имеет только одно значение в ячейке. Херт и Николз /12/ сформулировали следующее уравнение перемещения:
—С+£/*гя*С = 0 (2)
ы
где: С - объемная концентрация УОР скаляра, и V- вектор скорости.
-с I &• 1 р
Уравнение неразрывности:
^ + Vpi? = 0 dt У
(3)
может быть переписано в следующем виде:
-£-(logeP)+Vi? = 0
(4)
где р - плотность, и D/Dt полная производная по плотности. Допускаем, что плотность жидкости и воздуха постоянна, тогда уравнение (4) принимает вид:
где величина С может имеет только два значение: 1 - для жидкости и 0 - для воздуха.
Помимо описания поведения УОБ скаляра, для окончательного замыкания системы уравнений (1) необходимо определиться с математической моделью для описания турбулентных эффектов в потоке жидкости. В явном виде за их присутствие в системе уравнений (1) отвечает величина V. Согласно современным представлениям о природе турбулентных процессов, эффективная вязкость, является алгебраической суммой кинематической и турбулентной вязкости (последнюю еще называют коэффициентом турбулентного обмена). Причем если первая величина является константой жидкости, то по поводу теоретического определения второй величины, имеется большое количество мнений. Наиболее проверенным является так называемый к~Е подход, которым и воспользуемся. Согласно этому подходу, коэффициент турбулентной вязкости определяется как/9,13,7/: . .
где к - кинетическая энергия турбулентности; е - коэффициент диссипации.
Причем, в зависимости от режима потока, рекомендуется дополнительно пересматривать форму уравнений. Учитывая опыт стажировки автора в ВА\У (германском институте гидротехнического строительства), и результаты которые были получены в процессе совместной работы с немецкими коллегами над аналогичными задачами, в дальнейшем воспользуемся математической моделью для потоков с высокими числами Рейнольдса / 7 /. Уравнения для переноса турбулентной энергии к, тогда будет:
VU = О
Подставляя его в уравнение (2) получаем выражение:
—C + VCt/ = О
(5)
8t
(6)
(8)
где
(9)
g< 1 dp VbJ Р &i
(10)
Величина Рт- = 0, так как принимаем линейную зависимость между скоростями потока и турбулентными напряжениями в потоке. Скорость диссипации кинетической энергии е описывается по формуле:
1 д I г
«ре)+
А.
дх,
ри.е-
з Л
<7, дх,}
2 Г ди> ^ Эи,
]дх,
+с,
(П)
-с.гр~+с.
ре
ди, дх,
сте>Се1>^е2>^еЗ
где и -ы _ эмпирические коэффициенты, чьи значения, взятые из /
7 /. Следует обратить внимание на то, что традиционно поле турбулентности описывается через параметры * - Е. На практике же измеряют их две инварианты: / -интенсивность турбулентности (в процентах) и Ь - характеристическое расстояние (в метрах). Взаимосвязь между этими инвариантами пар величин осуществляется через следующие выражения:
к -1.5 * (V* I)
е = С„ 0.75 * к‘
71
(12)
(13)
где V - скорость потока;
Сц - константа модели (0.09).
Как правило, подходы в шлюз со стороны верхнего бьефа представляют собой небольшие участки с двух сторон огороженные бетонными палами. С одной стороны палы примыкают к верхней голове, а другой выходят в бьеф, где сечение потока увеличивается. Причем это увеличение может быть значительным (если имеем дело с так называемым озеровидным типом примыкания) или не значительным (если мы имеем дело с канальным типом примыкания). Часть типично канального примыкание к верхней голове представлена на рис. 1.
Рис. 1 Верхний, подходной канал с канальным типом примыкания
Для создания математической модели поведения свободной поверхности на базе уравнений (1) и (6), (8), (11), будем отталкиваться от области интегрирования КУРСОРЫ представленную на рис. 1. Помимо этого область интегрирования дополняется следующими граничными условиями:
1) сечение А ВЕР через которое поток попадает в камеру:
Ух = У(г) * соб (Ь) (14)
Уу = У(0 * ли (Ъ) (15)
К, = 0 (16)
¿ = 0.01 (17)
/=0.01 (18)
уог- 1 (19)
где У(/) - нормальная составляющая скорости потока в выходном сечении, которая определяется выражением:
У(0 = 0(0/8сеч (20)
где Q(t) - гидрограф наполнения камеры шлюза, являющийся гидравлической характеристикой процесса наполнения камеры шлюза.
8сеч - площадь сечения открытого отверстия, через которое поток попадает в камеру шлюза. Это сечение образуется между порогом шлюза и нижней кромкой рабочих ворот. Обе эти величины являются функциями от времени и зависят от эксплуатационных характеристик шлюза.
2) сечение {Ж/Р - область границы потока подходного примыкания к бьефу. Местоположение этого сечения назначается на расстоянии где уровень свободной поверхности воды будет колебаться с амплитудой близкой к нулю:
Р = Уводы* г + Рапш„ где 0<2<И, (21)
¿=0.01 (22)
/=0.01 (23)
УОР = 1 (24)
3) Сечение ¿МОС - горизонтальная плоскость ограничивающая расчетную область сверху:
К„ = 0 (25)
(26)
(27)
(28) (29)
4) для всех остальных граней расчетной области не протекания:
К = 0
к = о
е = 0
ёУОБ
-------= 0
с!п
Для задания начальных условий (в момент времени ¡—О) выделим две области заполняемые водой и воздухом соответственно:
1) для 0 <2<Н1
Р = уводы *2 + Ратм., где 0 < 2 < А; (34)
Ух=Уу=У2 = 0 (35)
¿ = 0.0 (36)
/=0.0 (37)
УОР = 1 (38)
— применим стандартное условие
(30)
(31)
(32)
(33)
¿ = 0.01 /=0.01 УОР = о
2) для HI < Z < Н2
Р - Ратм., где0<Z<ht (39)
1! II II О (40)
1 = 0.0 (41)
/=0.0 (42)
VOF-O (43)
Примечательной особенностью этой модели, является то, что с точки зрения математической физики данная модель будет относиться к классу задач с подвижными граничными условиями. Так как в процессе движения рабочих величина а (высота поднятия ворот над порогом) изменяется в интервале от 0 до 2.2 метра, то и граничная область ABEF не является постоянной. Причем, величина а является функцией от времени и индивидуальна для каждого шлюза.
Так как в качестве метода решения будет использоваться VOF метод контрольного объема, то для разрешения проблемы с подвижными граничными условиями, переформулируем задачу, как модель с граничными условиями первого рода. Такой путь вполне возможен, так как пространственное местоположение расчетной области неизменно, ибо меняются только место приложения и размеры граничной области ABEF и BCDE, сообразно процессу маневрирования затвором. Таким образом к граням тех ячеек, которые в данный конкретный момент времени принадлежат области выхода потока ABEF, применяются граничные условия (14) - (19) вместо условий не протекания (30) -(33). Такой прием позволяет максимально корректно смоделировать процессы волна формирования на данном участке без дополнительных приемов интегрирования.
На практике это будет выглядеть следующим образом. Имея график движения затвора - мы определяем величину открытия ворот а как функцию от времени, в любой момент времени. Помимо этого мы имеем величину расхода воды с гидрографа наполнения камеры и скорость течения воды в камеру V(t), для формул (14) и (15).
Для определения угла, под которым поток входит в камеру шлюза - Ь, необходимо иметь в виду, что он является функцией отношения глубины потока на пороге и величины открытия ворот а /I, 2, 10/. Но принимая во внимание результаты работы В.А. Кривошея можно с уверенностью утверждать, что траектория движения струи будет определяться напором только первые 10-20 сантиметров. Далее контуры потока попадающего в камеру шлюза будут определяться формой струенаправляющего козырька рабочих ворот. Подьемно-опускные ворота шлюзов ВДСК имеют уголковый струенаправляющий козырек, под углов в 45 градусов и высоту подъема ворот в 2,2 метра. Поэтому в дальнейших расчетах мы и будем оперировать именно этой величиной.
Естественно, что для численной реализации данной математической модели необходимо использовать соответствующее программное обеспечение и достаточный ресурс ПЭВМ. Как показали результаты стажировки автора в BAW, всем этим требованиях в полной мере удовлетворяют решатели систем COMET и STAR_CD /3, 5, 14, 15, 16/, которые предполагаются в дальнейшем к активному использованию для проведения исследований.
Список литературы
[1] Кривошей В.А. Гидродинамические нагрузки, действующие на подьемно-опускные ворота шлюзов, автореферат диссертации на соискание ученой степени к. т. н., Ленинград: 1985 г.
[2] Кривошей В.А. Исследование кинематики потока при истечении воды из-под ворот шлюза с головной системой питания // Сборник научных трудов. ЛИИВТ - гидротехнические сооружения и путевые работы на внутренних водных путях для судоходства. - Л.: ЛИИВТ, 1984. - С. 42-50
13J Методология по программе STAR CD. Computayional Dynamics Limited 1999.
[4] Патанкар С. Численные методы решения задач теплообмена и динамики жидкости / Пер. с англ. - М.: Энергоатомиздат, 1984. - 152 с.
[5] Руководство пользователя по программе STAR_CD. Computayional Dynamics Limited 1999
[6] Технический отчет по теме — “Bericht zu COMET numerische 3D-Simulation mit freier Wasseroberfläche” - Р.Лаузен - BAW - 99.3.03.06255.02 Карлсруэ. 2001 год.
[7] El Tahiy, S.H. 1983. 'k-e. equation for compressible reciprocating engine flows’, AIAA J. Energy, 7, No. 4, pp. 345-353.
[8] Gutachten über die seitliche Einleitung von Überschusswasser in den Einfahrtsbereich der Schleuse Nürnberg - технический отчет BAW. Автор Кастен Торенц - Карлсруэ - Bundesanstalt für Wasserbau (BAW) - Nr. 3.03.10043.00 - Dezember 2003.
[9] Launder, B.E., and Spalding, D.В. 1974. ‘The numerical computation of turbulent flow’, Comp.
Meth. in Appl. Mech. & Eng., 3, p. 269. ■
[10] Mises R. Berechnung Von Ausffluss und uberfallsahlen. V.D.I. № 22, 1917, p. 74-94.
[11] Harlow, F.H. and Welch, J.E. 1965. ‘Numerical calculation of time-dependent viscous incompressible flows of fluid with free surface’, Phys. Fluids, 8, pp 2182-2187.
[12] Hirt, C.W. and Nicholls, B.D. 1981. ‘Volume of Fluid (VOF) method for dynamical free boundaries’, J. Comput. Phys., 39, pp.201-225. >
[13] Rodi, W. 1979. ‘Influence ofbuoyancy and rotation on equations for turbulent length scalc’, Proc.
2nd Symp. on Turbulent Shear Flows.
[14] WWW.adapco.com
[15] WWW.cd.co.uk [16} WWW.starcd.uk
THE FORMULATION OF FULL 3D MODEL OF DYNAMICS DEVELOPMENT OF WAVE PROCESS IN THE TOP APPROACH CHANNEL OF A LOCK
I. Lipatov
The full three-dimensional mathematical model of development of wave process in the top approach channel of lock is formulated in this article. The model describes behavior of free surface and turbulent processes into flow. For model are formulated boundary and initialization conditions, model of turbulence. The end of article have the choice of the software and a necessary computing resource for practical realization of model.
УДК 502:531.001.57
Е. А. Лукина, к. т. н., ассистент.
Е. Ю. Чебан, аспирант, ВГАВТ.
603600, Нижний Новгород, ул. Нестерова, д. 5.
МОДЕЛИРОВАНИЕ ОБТЕКАНИЯ НЕФТЕСБОРНОГО БОНОВОГО ЗАГРАЖДЕНИЯ
Проблема уноса нефти под боновое заграждение может быть решена применением нефтесборного бонового заграяедения, оснащенного встроенным эжектором и приемным устройством. При определении характеристик этого устройства были численно решены уравнения Навъе-Стокса для потока жидкости при различных скоростях Представлены результаты моделирования обтекания бона при различных скоростях течения и расходах через приемное устройство.
Эффективность использования боновых заграждений для локализации разливов нефти зависит от множества параметров, среди которых основными являются скорость течения потока, вязкость нефтепродукта и осадка бона. В настоящее время раз-
