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

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

Том 18, № 4, 2013

Математическое моделирование эволюции берегового склона в каналах с песчаным руслом*

Б. В. Бондаренко, И. И. Потлпов Вычислительный центр ДВО РАН, Хабаровск, Россия e-mail: bvbondarenko@gmail.com, potapovii@rambler.ru

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

Ключевые слова: влекомые наносы, транспорт наносов, береговая эрозия.

Введение

Изучение процесса эволюции береговых склонов канала — береговой эрозии под действием протекающей в канале воды — является важной прикладной задачей. Первые попытки построения моделей береговой эрозии были сделаны в работах [1-5], в которых были предложены эвристические модели, позволившие получить качественные оценки изучаемого процесса. Дальнейшее развитие теории эрозии береговых склонов связано с большим количеством экспериментальных работ. Не претендуя на полноту, отметим некоторые из них [6-11]. Эти работы позволили сформулировать более совершенные теоретические модели, учитывающие механизмы транспорта влекомых, взвешенных и лавинных наносов при изучении процессов береговой эрозии [12-15]. Однако применение феноменологических формул транспорта донного материала позволило лишь качественно описывать исследуемые процессы. В работах [16-18] были предложены и развиты аналитические формулы транспорта донного материала, не содержащие феноменологических коэффициентов. Проверка работоспособности этих формул показала [19-21], что получаемые по ним расчётные данные хорошо согласуются с экспериментом [9-11].

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

*Работа выполнена при поддержке фонда ДВО РАН, грант X9 12-Ш-А-03-034, и РФФИ, грант № 12-01-98518-р восток(а).

риала определён в рамках модели, описанной в работе [16]. Лавинный механизм транспорта наносов, учитывающий обрушения берегового склона, определён в рамках модели [20]. Представлен алгоритм решения задачи, основанный на расщеплении задачи по физическим процессам. Выполнено сравнение результатов расчётов береговой эрозии для канала трапециевидной формы с известными экспериментальными данными [22].

1. Постановка задачи

Процесс деформации берегового склона рассмотрен на примере трапециевидного в начальный момент времени русла, имеющего постоянный малый продольный уклон 3 в направлении движения гидродинамического потока, движущегося в поле силы тяжести. В силу симметрии данного русла расчёт эволюции поперечного сечения канала выполнялся только для левой половины сечения. Геометрия расчётной области П представлена на рис. 1. Границы Г1 и Г4 подвижны и представляют собой смоченную береговую и донную поверхности и сухой береговой склон соответственно; Г2 — свободная поверхность потока, изменяющаяся во времени вследствие изменения геометрии донной поверхности; Г3 — граница, совпадающая с центром симметрии канала, изменяющаяся вследствие изменения глубины потока. Определение геометрии подвижных границ Г во времени — цель решения задачи.

Для вычисления границ Гд предложена математическая модель, содержащая уравнение движения установившегося гидродинамического потока

д2и д2и ^ - + — + Яд

дх\

дх2

0,

~ Рш9а „ где Яд =-, уравнение Экснера

№

ч д( до

(1 - с)р> ж + до = 0

и уравнение транспорта влекомых наносов [16]

о

Сот 3/2

( X А

tg р сое 7(1 - х) - - 2)

1 - Х> 0, 1 - X < 0,

0

Рис. 1. Геометрия и границы канала

C 4 p

Со

3 kyfpWFa tg * cos Y 9 k2dFa pro

то -' X = w

8 cx у т

Здесь x\, x2 — декартовы координаты расчётной области; u = u3(xi,x2) — скорость движения гидродинамического потока вдоль оси x3, направленной по потоку; — турбулентная вязкость потока, принимаемая равной постоянной величине; е — пористость данного материала; Fa = (ps — pw )g tg * cos y; т0 — напряжение начала движения донных частиц; * — угол внутреннего трения донных частиц; y — острый угол между нормалью к поверхности смеси и вертикальной линией; k = 0.4 — постоянная Кармана; cx — коэффициент лобового сопротивления частиц; d — диаметр частиц донного материала; pw, ps — плотность воды и донного материала; g — ускорение свободного падения; a — постоянный продольный уклон канала в направлении оси x3; t — однонаправленная временная координата; Z = Z(t, s) — кривая, определяющая профиль канала, включающий в себя дно русла и его недеформируемую береговую область; т — придонное касательное сдвиговое напряжение, определяемое в произвольной точке поверхности Z

как .

du

т = Цг

' dxi

Пг, (4)

С

где пг — компоненты вектора, нормального к донной поверхности £• Лавинный расход наносов оценивался по формуле [19, 22]

Ga = mps (1 — е)

D2 — tg2 *

m = ^ 11, Y > * 0, Y <

Однако прямо данная формула в расчётах не используется в силу сильного ухудшения жёсткости системы. Для учёта лавинного движения наносов использовалась модифицированная лавинная модель [23, 24], алгоритм которой приведен ниже. Задача, определяемая уравнениями (1)-(4), замыкается граничными

и(хг) = 0, хг Е Гь (6)

пгди = 0, хг Е Г2, Гз, (7)

С (*, 0) = Сс(0), (8)

(г,ь0)

ds

0

и начальными условиями

Z(0, s) = Zn (s), (10)

где ni — компоненты нормали к границе области П.

Изменение донной поверхности Z под воздействием потока приводит к изменению геометрии расчётной области П, что при постоянном руслоформирующем гидродинамическом расходе потока

Qn = / udn = const (11)

s

приводит к изменению гидродинамического сопротивления и уровня свободной поверхности потока Н.

Таким образом, решение эволюционной задачи развития профиля водопропускного канала заключается в решении гидродинамической задачи (1), определении придонных сдвиговых напряжений (4), по полученным значениям которых из решения уравнений (2), (3) находятся донные деформации. Далее с помощью дополнительного условия (11) следуют расчёт уровня свободной поверхности, соответствующего изменённой геометрии профиля канала, и циклическое повторение вышеуказанных шагов.

2. Метод решения гидродинамической задачи

и нахождения поверхностных сдвиговых напряжений

Для решения уравнения (1) и определения придонных сдвиговых напряжений т в работе использован прямой метод граничных элементов в постановке [25]. При этом граница Б области П разбивалась на N кусочно-постоянных элементов, а область П — на М треугольных симплекс-элементов, необходимых для вычисления функции и в каждой отдельной внутренней ячейке области и далее для вычисления расхода (11); д и р — индексы граничных элементов.

Отметим, что значительный временной вклад в решение задачи методом граничных элементов определяется затратами на расчёт правой части уравнения (1). С учётом того, что член Яь в уравнении (1) постоянен, последнее приводилось к однородному виду

д2и д2 и

дх\ + дх2

0

путём замены

и

и - Яд Х/2.

При этом граничные условия (6)-(7) преобразовывались к виду

и(жг) = Ядх2/2, Ж е Г1, ди

д%2 ди

0, ж; е Г2

дх1

= ЯдЖ1, Ж; е Гз.

'12) 13)

'14) :15) 16)

Согласно [25] для задачи (12), решаемой прямым методом граничных элементов, справедлива следующая система линейных уравнений:

ррч ич _ о™ Т1 = 0

:17)

где и1 и Т1 — векторы-столбцы размерности N, содержащие значения осевой скорости и гидродинамического потока и её нормальную производную Т в центрах граничных элементов соответственно; Гря,Орч — матрицы размерности N х N. Кодиагональные элементы матриц Ерч и Орч определяются как

Орч = — -

1

[г В (1п |г| - 1) + вк]

в В ,ГВ В

, ррч = _

вл,гл

2п

:18)

вл

га = V (жР - жА)2 + (ж2 - жА)2, гв = ^/(ж? - хВ)2 + (ж2 - жВ)2, (19)

к = ^ + + С|, (20) ^А2 + В2

где

-А — жр жв , В — ж^в жр, СС — жв (ж а ж в) ж в (жр жв ),

(жА, жА) и (жВ, жв) — координаты узлов д-го граничного элемента; (жр, жр) — координаты точки наблюдения, находящейся в геометрическом центре р-го граничного элемента,

7тр7Ар 7тр>Ар\ / 7тр7Вр 7тр7Вр

11 ^ -12 л , Вв = аг^ап ( 11 2 - ^ В

/тр/Ар + /тр/Ар I ' В I /тр/Вр + /тр/Вр

11 2 2/ \ 11 22

где

ВА = аг^ап -, вВ = аг^ап -, (21)

А I /тр/Ар + /тр/Ар / ' В — —

ж^ жр, /р жр жр, /р жр жр, ъ 1,2, (22)

= жАА'2 + ж?в2 - АВ (ж* - жА) = А (23)

ж1 = а2 + В2 , ж2 = ж2 в V 1 ж1 / . ^ '

Диагональные элементы матриц ^чр и Очр вычисляются по формулам

очч = -— (1п(Ь) - 1) , ^чч = 1, Ь = ^(жА - жВ)2 + (жА - жВ)2. (24) п^ 2 2

Решая систему линейных уравнений (17) методом Гаусса с выбором главного элемента, получаем полные векторы ич и Тч, определяя таким образом значения осевой скорости гидродинамического потока и и поверхностного сдвигового напряжения т на границе Б.

Скорость движения потока иг в г-м внутреннем узле области П вычисляется по формуле

иг = ич - ОгчТч - Яд(ж?)2/2. (25)

3. Метод решения задачи донных деформаций

Для решения задачи донных деформаций (2), (8)-(10) использован метод контрольных объёмов [19] (рис. 2), дискретный аналог которого для уравнения (2) имеет вид

«кС+1 = ад-иСЙ? + «^С-! + «к СП, (26)

где

_ Ак- _ Ак+ 0 _ А/к _ , ,0

«д-1 = ^—, «д+1 = ^—, = тг, «д = «д-1 + «д+1 + , (27) о/к- Аг

А/к = (ж1к+ - ж1к-)2 + (ж2к+ - ж2к-)2, (28)

= 0/(к+1)- = \! (ж1^+1 - ж1к )2 + (ж2к+! - ж2к )2. (29)

Ah

Рис. 2. Определение значений координат x2i в геометрических узлах

При выводе дискретного аналога полагалось, что донная поверхность £ = С(жьЖ2) в системе координат вк определена на каждом контрольном объёме как £ = С(в) и справедлива следующая линейная аппроксимация дна (н:

Zh

s \ s

1 - ÄT + ÄT Ö/k / Ö/k

k+i

или

где

1

cos Yk s

Xi

k+

x1k

X2k_ +

cos Yk s

-X2

X1k+ X1k-J

k+

cos Yk =

(X1k+ X1k-)

\l (X1k+ - X1k-)2 + (X2k+ - x2k- )

(30)

(31)

(32)

Формулы (30) или (31) позволяют определить начальные условия дискретной задачи (26)-(29) по донной поверхности £, заданной в системе координат (ж1,ж2).

После вычисления значений донных деформаций Д£к = Ск+1 — СП в каждом узле к-го контрольного объёма текущего п + 1-го временного слоя для определения значений координат х2к в геометрических узлах контрольных объёмов используется формула перехода

X2k

X2

k+

A(k

+

ACk+1 N

cos Yk+1 /

(33)

cos Yk

при этом горизонтальные координаты x1k остаются фиксированными в силу малости общих углов уклонов Yk •

При решении дискретной задачи (26)-(29) на каждом шаге по времени для каждого k-го контрольного объёма, для которого выполнялось условие Ak > 0, контролируется ограничение устойчивости схемы по времени

At <

min (A/k)

(34)

2 min (Ak)

Необходимо отметить, что в отличие от работы [19] расчёт донных деформаций и определение поля скоростей проводятся на смежных сетках. Это позволяет отказаться от мортар-сшивки сеток гидродинамической и русловой задач, ранее необходимой

h

Z

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

Для выполнения условия постоянства расхода потока (11) уровень свободной поверхности необходимо пересчитывать на каждом п-м шаге по времени. Коррекция уровня свободной поверхности выполняется в рамках итерационного процесса (к), описанного в [19]:

' 00,

= вгнп+Yr тг0 (нп - cm)+cm ), (35)

где

в = йг = АГ

г аг + АГ, аг + АГ,

НП+1, НП — уточняемые уровни свободной поверхности речного потока на п-м временном шаге; Як = / «*« — значение расх°да ре™ ^ нолученное для ур°вня

свободной поверхности НП; Я0 — заданный расход речного потока; П^ — площадь живого сечения потока на к-й итерации; АГ — временной шаг; (Ц = тт(£га) — минимальный уровень дна канала; аг — коэффициент релаксации схемы, принимающий значение, близкое к АГ.

Учёт лавинного движения наносов на береговом склоне производился по алгоритму, предложенному в работах [23, 24]. Для этого на каждом временном слое после вычисления значений (Ц выполнялась лавинная коррекция, при которой для каждого т-го узла контрольно-объёмной сетки проводились следующие парные вычисления:

ст = ст - ст-1 = ст-1 + (36)

где значение корректора Аь для левого берегового склона вычисляются как

Al Ч 0,+ас' Г> -Ас, г = Лж tg^ Acm = cm - G-i. (37)

Для предотвращения образования нефизичных осцилляций в основании берегового склона был применён метод нижней релаксации лавинной коррекции с параметром релаксации = 0.8.

4. Результаты численных расчётов

Для сформулированной модели, используя предложенный метод решения задачи, было выполнено численное моделирование эволюции берегового склона песчаного исходно трапециевидного канала. Расчёты проводились на адаптированной сетке, пример которой приведён на рис. 3, при следующих физико-механических и гранулометрических параметрах: к = 0.4, е = 0.4, = 1700 кг/м3, рад = 1000 кг/м3, д = 9.81 м/с2, 3 = -0.02, сх = 0.4, ^ = 0.002 м. Построение сетки выполнялось ОшзЬ-генератором [26] с адаптацией по диаметру элемента в направлении дна русла в среднем на один порядок.

Результаты численных расчётов и данные натурного эксперимента [22] представлены на рис. 4. Линии 2- 4 получены в результате расчёта по приведённому алгоритму. Линия 1 определяет начальную геометрию канала. В начальный момент времени ¿0 = 0 согласно [22] в канале задан гидравлический расход потока Q = 13 л/с. В моменты времени ^ = 250 и ¿2 = 760 мин Q увеличивается соответственно в 1.5 и 2.5 раза относительно начального расхода. При этом свободная поверхность потока поднимается и происходит затопление береговой части русла. Экспериментальные данные [22] в различные моменты времени (240, 570, 820 мин) представлены точками.

Из приведённых графиков 2- 4 видно, что береговой склон в различные моменты времени имеет чётко выраженный лавинный участок берега. Уклон данного участка близок к углу естественного откоса донного материала Отметим: если не учитывать предложенных в алгоритме лавинных механизмов движения наносов, то получаемые решения теряют физический смысл — движение берега прекращается, а его уклон в точке начала движения донного материала стремится к 90 °. В процессе отступления берегового склона лавинное движение донного материала в своей глубоководной части переходит в движение, обусловленное турбулентно-диффузионным механизмом переноса донного материала (определяемого уравнением (3)). При этом происходит замывание глубоководной части канала. Расширение глубоководной части дна происходит с увеличением его уклона в сторону берега и уменьшением лавинной части берегового склона. Изменение донного профиля и общий подъём уровня дна приводят к небольшому (менее 5 % от исходного) подъёму уровня свободной поверхности потока. Основной подъём уровня свободной поверхности происходит при увеличении гидравлического расхода потока в моменты времени ¿2.

Рис. 3. Пример адаптивной расчётной сетки, полученной при решении гидродинамической задачи

Г,

Ах2

0.14

• ? = 240 мин

0.12

0.1

0.08

0.4

0.6

0.8

Рис. 4. Эволюция донной поверхности ( во времени

Численное решение, полученное в рамках предложенной модели, показало хорошее согласование расчётных и экспериментальных профилей — относительная погрешность расчётных данных в среднем не превышает 10%, что сопоставимо с системной точностью выполненных экспериментов. Расчётная и экспериментальная скорости отступления берегового склона различаются не более чем на 6 %. Максимальное различие этих данных, достигающее 47%, наблюдается на поздних временах счёта в месте сопряжения лавинной и гравитационной областей движения наносов, после увеличения гидравлического расхода в момент времени t2. К сожалению, для определения причины такого отклонения имеющихся экспериментальных данных недостаточно. Полученное различие расчётных и экспериментальных данных может быть обусловлено следующими причинами (или их совокупностью):

— нестационарностью потока в периоды изменения гидравлического расхода;

— выбранной турбулентной моделью потока;

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

Заключение

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

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

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

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

[1] Glover R.E., Florey Q.L. Stable Channel Profiles. U.S. Bureau of Reclamation. Washington, 1951.

[2] Макавеев Н.И. Русло реки и эрозия в её бассейне. М.: Изд-во АН СССР, 1955. 348 с.

[3] Русловой процесс / Н.Е. Кондратьев, А.Н. Ляпин, И.В. Попов, С.И. Пиньковский, Н.Н. Федоров, И.И. Якунин. Л.: Гидрометеоиздат, 1959. 372 с.

[4] Callander R.A. Instability and river channels // J. Fluid Mech. 1969. Vol. 36, No. 3. P. 465-480.

[5] Ikeda S., Parker G., Saway K. Bend theory of river meanders. Part 1. Linear development // Ibid. 1981. Vol. 112. P. 363-377.

[6] Muneo Hirano. River-bed variation with bank erosion // J. of Hydraulic, Coastal and Environmental Eng. (in Japanese). 1973. No. 210. P. 13-20.

[7] Shoji Fukuoka, Masashige Yamasaka. Shear stress distribution on a continuous boundary and widening process of a straight alluvial channel // Japan Soc. of Civil Eng. 1984. No. 351, II2. P. 87-96.

[8] Mtivier F., Meunier P. Input and output mass flux correlations in an experimental braided stream. Implications on the dynamics of bed load transport // J. of Hydrology. 2003. No. 271. P. 22-38.

[9] Пушкарёв В.Ф. Движение влекомых наносов // Труды ГГИ. 1948. Вып. 8 (62).

[10] Шуляк Б.А. Физика волн на поверхности сыпучей среды и жидкости. М.: Наука, 1971.

[11] Engelund F., Fredsoe J. Sediment ripples and dunes // Annual Rev. of Fluid Mech. 1982. No. 14. P. 13-37.

[12] Parker G. Self-formed straight rivers with equilibrium banks and mobile bed. Part 1. The sand-silt river //J. Fluid Mech. 1978. Vol. 89. P. 109-126.

[13] Ikeda S. Stable channel cross-sections of straight sand rivers //J. Water Resources Res. 1991. Vol. 27, No. 9. P. 2429-2438.

[14] Diplas P., Vigilar G. Hydraulic geometry of threshold channels //J. Hydr. Eng. ASCE. 1992. Vol. 118, No. 4. P. 597-614.

[15] Farias H.D. Physical and mathematical modelling of alluvial channels in regime, HYDRA 2000 // Hydraulics Res. and its Appl. Next Century / Eds. A. Ervine, T. Telford. London, U.K., 1995. Vol. 1. P. 348-353.

[16] Петров П.Г. Движение сыпучей среды в придонном слое жидкости // ПМТФ. 1991. Т. 41, № 5. С. 72-75.

[17] Петров А.Г., Петров П.Г. Вектор расхода наносов в турбулентном потоке над размываемым дном // Там же. 2000. Т. 41. С. 102-112.

[18] Петров А.Г., Потапов И.И. О развитии возмущений песчаного дна канала // Докл. АН. 2010. Т. 431, № 2. С. 191-195.

[19] Потапов И.И., Бондаренко Б.В. Моделирование эволюции поперечного сечения песчаного канала // Вычисл. технологии. 2009. Т. 14, № 5. С. 1-14.

[20] Потапов И.И., Щекачёва М.А. Определение скорости размыва берегового склона в реках с песчаным дном // Вестник Удмуртского гос. ун-та. Механика. 2011. Вып. 4. С. 116-120.

[21] Петров А.Г., Потапов И.И. О механизмах развития донных волн в канале с песчаным дном // ПМТФ. 2011. Т. 52, № 2. С. 81-91.

[22] Pitlick J., Pizzuto J.E., Marr J. Width adjustment in alluvial channels. St. Anthony Falls Laboratory (http://www.colorado.edu/geography/geomorph/nsf_safl.html)

[23] Tjerry S. Morphological Calculations of Dunes in Alluvial Rivers. Ph.D. Thesis. ISVA, Techn. Univ. of Denmark, 1995. 193 p.

[24] Sanne L.N. Modelling of Sand Dunes in Steady and Tidal Flow. Ph.D. Thesis. MEK-DTU. Coastal ancl River Eng. Section. Tech. Univ. of Denmark, 2003. 185 p.

[25] Бенерджи П., Баттерфилд Р. Метод граничных элементов в прикладных науках: Пер. с англ. М.: Мир, 1984. 494 с.

[26] Geuzaine C., Remacle J.-F. Gmsh: A three-dimensional finite element mesh generator with built-in pre- and post-processing facilities // Intern. J. for Numer. Methods in Eng. 2009. Vol. 79, No. 11. P. 1309-1331.

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