Краевые условия и расчёт поверхностных напряжений в численных моделях естественных водотоков Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 18, № 5, 2013

Краевые условия и расчёт поверхностных напряжений

>к

в численных моделях естественных водотоков*

В. А. Шлычков Институт водных и экологических проблем СО РАН (Новосибирский филиал), Россия e-mail: slav@iwep.nsc.ru

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

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

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

1. Условия на дне

При достаточно большой скорости потока донные отложения, представляющие собой смесь твёрдых частиц разной крупности, становятся подвижными, что даёт начало процессу донных деформаций [1]. Условия предельного равновесия частиц формулируются в терминах касательных напряжений на дне т, причём критерий потери устойчивости записывается в виде т > тсг, где тсг — некоторое критическое значение сдвигового напряжения. Отсюда следует, что при решении задач перемещения донных наносов необходимо по возможности точно рассчитать величину донных напряжений в водной среде. Определение касательных напряжений крайне важно для расчёта динамики наносов в водотоках при прогнозах русловых деформаций и обеспечении условий судоходства, бесперебойной работы речных водозаборов и других гидротехнических сооружений.

* Работа выполнена в рамках Интеграционного проекта № 109 СО РАН.

В рамках гипотезы Буссинеска параметризацию напряжений Рейнольдса п'и)' можно записать в виде

-п^ = т = Кдп = п2, (1)

р д х

где К — коэффициент турбулентного обмена, п* — динамическая скорость [2], характеризующая масштаб турбулентных пульсаций, х — вертикальная координата. Экспериментальное и теоретическое определение коэффициента К представляет собой одну из основных задач теории турбулентности.

Введённый Прандтлем "путь перемешивания" со временем утратил свое значение, и современные теории используют другие гипотезы замыкания. Так, в [3] для турбулентного ядра в пристеночном пограничном слое предложена зависимость

К = кп*х,

которая приводит к логарифмическому распределению скорости по глубине

п 1 , х

— = -1п—+ Б3, х > хо, (2)

п* к х0

где к ~ 0.4 — постоянная Кармана, Б8 — постоянная интегрирования, подлежащая определению. Экспериментально установлено, что Б8 не является постоянной, а представляет собой функцию (пристеночная Б-функция [4]), зависящую от режима обтекания поверхности. Возникает задача построения Б-функции для реальных водотоков. Асимптотически с возрастанием локального числа Рейнольдса Б8 стремится к постоянному значению 8.48, характеризующему отношение п0/п* в опытах Никурадзе (здесь п0 — скорость потока на уровне выступов шероховатости х0). Асимптотическое значение Б8 часто применяется в расчётах [5], однако является довольно грубым и не согласуется с подходами, принятыми в классической гидравлике речных русел. Чтобы показать это, рассмотрим уравнение движения для равномерного потока, однородного в продольном и поперечном направлениях, с ровным шероховатым дном и уклоном гх

дп . д т^дп

дь = + д~хКд~х. (3)

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

дп дп

К—— = сйп|п| при х = х0, К—— = 0 при х = (, (4)

д х д х

где величина сопротивления Сй задаётся безразмерной, х = £ — уровень свободной поверхности. Физически корректное задание параметра с^ является важным условием адекватного воспроизведения профиля скорости и, согласно (1), донных касательных напряжений. С учётом логарифмичности (в виде профиля скорости, подчиняющегося зависимости (2)) из первого условия в (4) получаем

1

Сй = Б. (5)

8

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

С целью вывода гидравлически приемлемого выражения для са проинтегрируем уравнение (3) по вертикали от г0 до £ и обозначим через К = £ — глубину потока. С учётом краевых условий (4) получим

^ = — дК1х — сап о|ио|, (6)

С

где q = ! идьХ — удельный расход воды.

¿0

В гидравлических исследованиях речных потоков осреднённое по глубине уравнение движения традиционно имеет вид [1]

dt = -9hix - u|u1, (7)

где U — средняя по вертикали скорость й = q/h, Cs — коэффициент Шези, характеризующий шероховатость русла. Коэффициент Шези табулирован для русел различных типов и задаётся исходя из общих характеристик водотока. Гидравлический аналог динамической скорости вычисляется в соответствии с (7) как U = дй2/С2. Таким образом, имеется двойственное представление для динамической скорости — следующее из модели с вертикальной детализацией (1), (4) и вытекающее из (7). В обоих случаях накоплен обширный теоретический и экспериментальный материал для идентификации этих параметров по отдельности, но сами подходы существуют порознь и не являются согласованными. Сопоставительный анализ их общей базы позволяет найти непротиворечивое определение коэффициента сопротивления и даёт конструктивный метод его расчёта.

Сравнивая уравнения (6), (7), получим условия их совместности в виде

д и |u|

cd = ~k~ гЛ. (8)

С2 Uo |uo|

Отношение U/u0 заранее не известно, так как зависит от искомого решения Z

- 1 f

(и = — udz). При этом характер нелинейности в первом условии (4) усложняется

zo

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

Значение U может быть вычислено также с помощью любой формулы, аппроксимирующей профиль скорости. Например, в предположении, что u(z) явно задано в виде (2), после аналитического вычисления U и его подстановки в (8) получим (при этом

малой величиной — « 10-4 пренебрегалось) h

cd

1 'ln - - 1|+1

кБА zo

2

L 2 .

s

C

Сс1

0.01 ■

0.02-

0 2 4 6 8 10 г0,мм

Рис. 1. Зависимость (8) коэффициента сопротивления са от диаметра наносов го при С = 50, 60, 70 (кривые 1-3); 4 — значение са, вычисленное по (5)

Формула (9) связывает локальную характеристику донной шероховатости (параметр г0) и макросопротивление русла (параметр С8, отражающий наличие на дне гряд, растительности, замусоренности ложа, эффектов меандрирования и др.). Вычисленное по (9) значение са для естественных русел находится в диапазоне 0.01-0.05, тогда как из представления (5) следует са ~ 0.014 (рис. 1, кривая 4).

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

Адекватный расчёт сверхвысоких градиентов логарифмической функции в придонной области требует сеточного разрешения, сопоставимого с вертикальным размером выступов шероховатости, т. е. долей миллиметра. В многомерных моделях это приводит к чрезмерно большому числу переменных, что не экономично с вычислительной точки зрения. По этой причине в численных моделях граничное условие на уровне шероховатости г0 обычно сносят на вышележащий счётный уровень , определяемый шагами сетки, который, тем не менее, должен располагаться не выше слоя логарифмичности (толщиной примерно 1/6 от глубины потока). С изменением уровня сопряжения изменится и значение коэффициента сопротивления, причём его новое значение cs может быть получено в виде

2. Условия на свободной поверхности

Современные подходы к описанию вертикального турбулентного обмена в водоёмах основаны на применении двухпараметрической модели для кинетической энергии турбулентности (КЭТ) е и скорости её диссипации е, которые в одномерном приближении имеют вид [6]

са

(10)

де д Ве — = ае—К— + КЗ - е, дЬ дг дг

де д т^де е ^ т е2 ...

— = К— + С22-КЗ - с3 —, (11)

дЬ дг дг е е

е2

К = с,-,

е

ди\2 ддр

где ае,ае,с,,с2,с3 — эмпирические постоянные, З = —--—--скорость попол-

дг р дг

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

и2 е3/2

"" * —Се —, (12)

V с,' г

где се — комбинация постоянных. Заметим, что хотя решение (12) справедливо лишь в придонной части потока, где выполняется гипотеза квазистационарного подслоя [2] де/дг = 0, оно может быть использовано как прообраз для формулирования условий в общей постановке.

Для системы (11) поставим краевую задачу. Граничные условия на дне зададим в виде соотношений, вытекающих из (12):

де е 2

— = 0, е = се— при г = го. (13) дг г0

На свободной поверхности при отсутствии ветра нередко используют условия [7]

де де

_ = °, _ = о ПрИ, = с (14)

Эти условия не обеспечивают затухания К с приближением к свободной поверхности, как это должно быть согласно экспериментальным данным [1]. Подавление турбулентного обмена здесь объясняется тем, что свободная поверхность в стационарном течении неподвижна в пространстве и её можно рассматривать как непроницаемую твёрдую стенку, на которой турбулентные пульсации затухают. Однако из системы (11) при условиях (14) вытекает д К/д г = 0, что означает наличие экстремума, т. е. максимальную интенсивность турбулентного обмена на свободной поверхности. Отсюда следует, что условия (14) физически не корректны и следует использовать другие условия, например [8]

де е3/2

— = 0, е = се- при г = С, (15)

дг га

где га — масштаб турбулентности на свободной поверхности.

На рис. 2 приведены профили и, К, е в потоке глубиной 5 м, полученные интегрированием уравнений (3), (11) на основе 150-уровенной (Дг = 3.3 см) численной модели методом установления. Высота выступов шероховатости г0 принята равной 0.2 мм. Другие параметры подбирались так, чтобы максимум скорости был равен 1м/с. Видно,

Рис. 2. Расчётные профили продольной скорости течения (1), коэффициента турбулентности (2) и кинетической энергии турбулентности (3)

что максимальные значения сдвига скорости du/dz наблюдаются на уровне шероховатости. Имея в виду оценку точности воспроизведения динамической скорости, укажем, что численное значение агрегата Kuz = ul оказалось равным 8.8см2/с2, а точное равно 11см2/с2, т.е. относительная погрешность составила около 20%. Прямое вычисление напряжения трения на дне связано с произведением очень малой величины K на очень большую величину du/dz, что и обусловливает генерацию ошибок при машинных расчётах.

Воздействие ветра в условиях открытой воды задаётся поверхностными ветровыми напряжениями та

TSdu та /■

к — = — при z = Z. (16)

dz р

При формулировании условия для КЭТ иногда используют зависимость [7]

/ \ 3/2

aeK— = cJ при z = (, (17)

где ca — постоянная. Формула (17) правильно отражает размерность левой части, но неверна с физической точки зрения. Чтобы показать это, рассмотрим уравнения энергии, полученные интегрированием по z уравнения движения (3), умноженного на u, и уравнение для КЭТ в (11). После подстановки краевых условий (полагаем u > 0, р = const) получим

дЕ _ / • \ /V fdu\

та

dt к{ж) / + 7ulz=z- Cd|uo1

3

д(е) ггде. /гг/'ди\2\ . . , ,

ж = ^^ + \к{ж) /-{е}- (18)

где угловые скобки означают интегрирование по вертикали, Е = и2/2 — кинетическая энергия среднего движения.

2

Рис. 3. Схема превращений энергии в модели

При сложении уравнений в (18) слагаемые, описывающие генерацию КЭТ, взаимно уничтожаются, а оставшиеся слагаемые имеют ясный физический смысл. Схема энергетических преобразований, показанная на рис. 3, включает источники, связанные с ветровыми напряжениями и уклоном дна, и стоки энергии в вязкий подслой и в кинетическую энергию турбулентности. Последняя диссипирует во всей толще потока. Если же учесть условие (17), согласно которому КЭТ-составляющая генерируется непосредственно из ветрового источника, то общий баланс будет нарушен вследствие формирования фиктивной скорости пополнения суммарной энергии системы (перечёркнутая связь на рис. 3). Если использовать предположение о прямом пополнении КЭТ за счёт обрушения волн, то ветровой источник следует разбить на два компонента, один из которых будет обусловливать формирование дрейфовой составляющей энергии, а другой — КЭТ. Этот подход используется, например, в [9]. Таким образом, и при наличии ветрового воздействия условия (15) остаются физически корректными.

Список литературы

[1] Гришлнин К.В. Основы динамики русловых потоков. М.: Транспорт, 1990. 320 с.

[2] Монин А.С., Яглом А.М. Статистическая гидромеханика. Т. 1. СПб: Гидрометеоиздат, 1992. 694 с.

[3] Ландау Л.Д., Лифшиц Е.М. Механика сплошных сред. М.: Гостехтеориздат, 1953. 736 с.

[4] Зеленецкий В.А., Богатко Т.В. Об одном способе описания пристенного слоя // Прикладная механика и техническая физика. 2011. Т. 52, № 2. С. 73-80.

[5] Игнатова Г.Ш., Квон В.И. О гидродинамической схеме скольжения при турбулентном течении // Метеорология и гидрология. 1978. № 7. С. 50-54.

[6] Роди В. Модели турбулентности окружающей среды // Методы расчёта турбулентных течений. М.: Мир, 1984. C. 227-322.

[7] Alfrink B.J., van Rijn L.C. Two-equation turbulence model for flow in trenches // J. of Hydraulic Eng. 1983. Vol. 109, No. 7. P. 941-958.

[8] Naot D., Rodi W. Calculation of secondary currents in channel flow // J. of Hydraulics Division. 1982. Vol. 108, No. 8. P. 948-968.

[9] Китайгородский С.А. Физика взаимодействия атмосферы и океана. Л.: Гидрометеоиздат, 1970. 290 с.

Поступила в 'редакцию 27 декабря 2012 г., с доработки — 8 апреля 2013 г.