Чисельне моделювання еволюції нерівностей великого масштабу на річковому дні
Виконане чисельне моделювання переносу домішок у придонній області відкритих водойм за наявності там нерівностей великого масштабу, яке грунтується на системі рівнянь мілкої води і рівнянні еволюції дна Екснера. Розвинений чисельний алгоритм поєднує центральну схему другого порядку Курганова-Ноелля-...
Gespeichert in:
| Datum: | 2015 |
|---|---|
| 1. Verfasser: | |
| Format: | Artikel |
| Sprache: | Ukrainian |
| Veröffentlicht: |
Інститут гідромеханіки НАН України
2015
|
| Schriftenreihe: | Прикладна гідромеханіка |
| Schlagworte: | |
| Online Zugang: | http://dspace.nbuv.gov.ua/handle/123456789/116497 |
| Tags: |
Tag hinzufügen
Keine Tags, Fügen Sie den ersten Tag hinzu!
|
| Назва журналу: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Zitieren: | Чисельне моделювання еволюції нерівностей великого масштабу на річковому дні / І.М. Горбань // Прикладна гідромеханіка. — 2015. — Т. 17, № 1. — С. 21-36. — Бібліогр.: 32 назв. — укр. |
Institution
Digital Library of Periodicals of National Academy of Sciences of Ukraine| id |
irk-123456789-116497 |
|---|---|
| record_format |
dspace |
| spelling |
irk-123456789-1164972017-04-29T03:02:58Z Чисельне моделювання еволюції нерівностей великого масштабу на річковому дні Горбань, І.М. Науковi статтi Виконане чисельне моделювання переносу домішок у придонній області відкритих водойм за наявності там нерівностей великого масштабу, яке грунтується на системі рівнянь мілкої води і рівнянні еволюції дна Екснера. Розвинений чисельний алгоритм поєднує центральну схему другого порядку Курганова-Ноелля-Петрової для інтегрування гідродинамічних рівнянь зі зваженою істотно неосцилюючою схемою "проти течії" п'ятого порядку, яка моделює транспорт наносів. Показано, що зв'язане моделювання гідро- та морфодинамічних процесів разом із застосуванням схем високого порядку забезпечує стійкість результатів у довгострокових розрахунках. Одержані часові і просторові характеристики еволюції піщаних кар'єрів та пагорбів на річковому дні, які показують, що в процесі розмиву не лише змінюється поперечний профіль нерівності, а й відбувається її перенос за течією. Выполнено численное моделирование переноса примесей в придонной области открытых водоемов при наличии там крупномасштабных неровностей. Оно основывается на системе уравнений мелкой воды и уравнении эволюции дна Экснера. Развитый численный алгоритм объединяет центральную схему второго порядка Курганова-Ноэлля-Петровой для интегрирования гидродинамических уравнений со взвешенной существенно неосциллирующей схемой "против течения" пятого порядка, которая моделирует перенос взвешенных частиц. Показано, что связанное моделирование гидро- и морфодинамических процессов, а также применение схем высокого порядка обеспечивает устойчивость результатов в долгосрочных расчетах. Получены временные и пространственные характеристики песчаных карьеров и холмов, которые показывают, что в процессе размыва не только изменяется поперечный профиль неровности, а и происходит ее перенос по течению. This paper deals with numerical simulation of bed-load sediment transport near large-scale irregularities placed on the bottom of open channels. It is founded on a shallow water system and a so called Exner equation that describes the evolution of topography. The numerical algorithm developed is coupled the Kurganov-Noelle-Petrova semidiscrete central-upwind scheme for integrating the hydrodynamic equations with a fifth order Euler-WENO scheme for sediment transport modeling. The splitting simulation of hydro- and morphodynamical processes as well as application of the high order schemes were shown to ensure stability of the results during a long-term modeling. The temporal and spatial characteristics of evolution of large-scale bottom irregularities such as sand humps and cavities were obtained. Those point out that sediment transport not only changes the cross profile of the irregularities but causes their moving on down the stream. 2015 Article Чисельне моделювання еволюції нерівностей великого масштабу на річковому дні / І.М. Горбань // Прикладна гідромеханіка. — 2015. — Т. 17, № 1. — С. 21-36. — Бібліогр.: 32 назв. — укр. 1561-9087 http://dspace.nbuv.gov.ua/handle/123456789/116497 532.517 uk Прикладна гідромеханіка Інститут гідромеханіки НАН України |
| institution |
Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| collection |
DSpace DC |
| language |
Ukrainian |
| topic |
Науковi статтi Науковi статтi |
| spellingShingle |
Науковi статтi Науковi статтi Горбань, І.М. Чисельне моделювання еволюції нерівностей великого масштабу на річковому дні Прикладна гідромеханіка |
| description |
Виконане чисельне моделювання переносу домішок у придонній області відкритих водойм за наявності там нерівностей великого масштабу, яке грунтується на системі рівнянь мілкої води і рівнянні еволюції дна Екснера. Розвинений чисельний алгоритм поєднує центральну схему другого порядку Курганова-Ноелля-Петрової для інтегрування гідродинамічних рівнянь зі зваженою істотно неосцилюючою схемою "проти течії" п'ятого порядку, яка моделює транспорт наносів. Показано, що зв'язане моделювання гідро- та морфодинамічних процесів разом із застосуванням схем високого порядку забезпечує стійкість результатів у довгострокових розрахунках. Одержані часові і просторові характеристики еволюції піщаних кар'єрів та пагорбів на річковому дні, які показують, що в процесі розмиву не лише змінюється поперечний профіль нерівності, а й відбувається її перенос за течією. |
| format |
Article |
| author |
Горбань, І.М. |
| author_facet |
Горбань, І.М. |
| author_sort |
Горбань, І.М. |
| title |
Чисельне моделювання еволюції нерівностей великого масштабу на річковому дні |
| title_short |
Чисельне моделювання еволюції нерівностей великого масштабу на річковому дні |
| title_full |
Чисельне моделювання еволюції нерівностей великого масштабу на річковому дні |
| title_fullStr |
Чисельне моделювання еволюції нерівностей великого масштабу на річковому дні |
| title_full_unstemmed |
Чисельне моделювання еволюції нерівностей великого масштабу на річковому дні |
| title_sort |
чисельне моделювання еволюції нерівностей великого масштабу на річковому дні |
| publisher |
Інститут гідромеханіки НАН України |
| publishDate |
2015 |
| topic_facet |
Науковi статтi |
| url |
http://dspace.nbuv.gov.ua/handle/123456789/116497 |
| citation_txt |
Чисельне моделювання еволюції нерівностей великого масштабу на річковому дні / І.М. Горбань // Прикладна гідромеханіка. — 2015. — Т. 17, № 1. — С. 21-36. — Бібліогр.: 32 назв. — укр. |
| series |
Прикладна гідромеханіка |
| work_keys_str_mv |
AT gorbanʹím čiselʹnemodelûvannâevolûcíínerívnostejvelikogomasštabunaríčkovomudní |
| first_indexed |
2025-07-08T10:29:23Z |
| last_indexed |
2025-07-08T10:29:23Z |
| _version_ |
1837074287989620736 |
| fulltext |
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
УДК 532.517
ЧИСЕЛЬНЕ МОДЕЛЮВАННЯ ЕВОЛЮЦIЇ НЕРIВНОСТЕЙ
ВЕЛИКОГО МАСШТАБУ НА РIЧКОВОМУ ДНI
I. М. Г О РБ А Н Ь
Iнститут гiдромеханiки НАН України, Київ
вул. Желябова, 8/4, 03680, МСП, Київ-180, Україна
E-mail: ivgorban@gmail.com
Одержано 30.01.2015
Виконане чисельне моделювання переносу домiшок у придоннiй областi вiдкритих водойм за наявностi там
нерiвностей великого масштабу, яке ґрунтується на системi рiвнянь мiлкої води i рiвняннi еволюцiї дна Екснера.
Розвинений чисельний алгоритм поєднує центральну схему другого порядку Курганова-Ноелля-Петрової для
iнтегрування гiдродинамiчних рiвнянь зi зваженою iстотно неосцiлюючою схемою "проти течiї"п’ятого порядку,
яка моделює транспорт наносiв. Показано, що зв’язане моделювання гiдро- та морфодинамiчних процесiв разом iз
застосуванням схем високого порядку забезпечує стiйкiсть результатiв у довгострокових розрахунках. Одержанi
часовi i просторовi характеристики еволюцiї пiщаних кар’єрiв та пагорбiв на рiчковому днi, якi показують, що в
процесi розмиву не лише змiнюється поперечний профiль нерiвностi, а й вiдбувається її перенос за течiєю.
КЛЮЧОВI СЛОВА: донна поверхня, розмив, перенос домiшок, заглиблення, пiщаний пагорб
Выполнено численное моделирование переноса примесей в придонной области открытых водоёмов при наличии
там крупномасштабных неровностей. Оно основывается на системе уравнений мелкой воды и уравнении эволюции
дна Экснера. Развитый численный алгоритм объединяет центральную схему второго порядка Курганова-Ноэлля-
Петровой для интегрирования гидродинамических уравнений со взвешенной существенно неосциллирующей схемой
"против течения"пятого порядка, которая моделирует перенос взвешенных частиц. Показано, что связанное мо-
делирование гидро- и морфодинамических процессов, а также применение схем высокого порядка обеспечивает
устойчивость результатов в долгосрочных расчётах. Получены временные и пространственные характеристики пе-
счаных карьеров и холмов, которые показывают, что в процессе размыва не только изменяется поперечный профиль
неровности, а и происходит её перенос за течением.
КЛЮЧЕВЫЕ СЛОВА: донная поверхность, размыв, перенос примесей, углубления, песчаный холм
This paper deals with numerical simulation of bed-load sediment transport near large-scale irregularities placed on the
bottom of open channels. It is founded on a shallow water system and a so called Exner equation that describes the
evolution of topography. The numerical algorithm developed is coupled the Kurganov-Noelle-Petrova semidiscrete central-
upwind scheme for integrating the hydrodynamic equations with a fifth order Euler-WENO scheme for sediment transport
modeling. The splitting simulation of hydro- and morphodynamical processes as well as application of the high order
schemes were shown to ensure stability of the results during a long-term modeling. The temporal and spatial characteristics
of evolution of large-scale bottom irregularities such as sand humps and cavities were obtained. Those point out that
sediment transport not only changes the cross profile of the irregularities but causes their moving on down the stream.
КЛЮЧЕВЫЕ СЛОВА: river bottom, erosion, sediment transport, hollow, sand dune
ВСТУП
Рiчкова течiя, як динамiчна система, складає-
ться iз взаємопов’язаних гiдрологiчного i морфо-
логiчного процесiв. Гiдродинамiчнi явища визна-
чаються глибиною i швидкiстю водного потоку, на-
явнiстю вихорiв i активних турбулентних зон то-
що. Морфологiчнi процеси включають розмив рi-
чкового дна i берегiв, втягнення утворених части-
нок у рух, їхнiй перенос i випадання осадiв. Про-
гнозування гiдродинамiчних i морфологiчних змiн
у вiдкритих водоймах є важливим для отримання
характеристик навколишнього середовища та йо-
го еволюцiї з метою попередження негативних на-
слiдкiв техногенних i природних катастроф.
Протягом довгого часу наукову базу для дослiд-
ження рiчкових течiй забезпечувала гiдравлiка. Її
концептуальнi моделi ґрунтуються на емпiричних
спiввiдношеннях, одержаних або зi спостережень
за полем течiї, або з лабораторних експериментiв
[1, 2]. В 90-х роках минулого сторiччя в рiчко-
ву гiдродинамiку почали широко впроваджувати-
ся комп’ютернi технологiї. Їхньою перевагою є те,
що вони реалiзують пряме масштабне моделюван-
ня течiї, дозволяють виявити вплив окремих фа-
кторiв i є набагато дешевшими за натурнi чи мо-
дельнi експерименти [3]. Крiм того, вони ґрунтую-
ться на строгих математичних моделях, що усуває
проблеми, пов’язанi з вибором великої кiлькостi
емпiричних залежностей.
Зазвичай, для чисельного моделювання неста-
цiонарних руслових потокiв з переносом домiшок у
придоннiй областi використовуються система рiв-
нянь мiлкої води, що описує водний потiк, i рiвня-
ння збереження маси наносiв, за яким розраховує-
ться еволюцiя дна. Складнiсть проблеми полягає в
тому, що часовi масштаби змiн у рiчковому потоцi
c© I. М. Горбань, 2015 21
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
i на днi суттєво вiдрiзняються, iнколи – на декiль-
ка порядкiв. Їхнє спiввiдношення необхiдно вра-
ховувати при побудовi чисельних алгоритмiв для
розв’язання математичних рiвнянь [4].
Якщо водний потiк близький до стацiонарного i
слабо залежить вiд змiн на доннiй поверхнi, гiдро-
морфодинамiчна система може розв’язуватися по-
слiдовно, коли спочатку iнтегруються гiдродина-
мiчнi рiвняння до досягнення течiєю стану рiвно-
ваги, а потiм оновлюється форма дна. Такий спро-
щений пiдхiд застосовується для розв’язання бага-
тьох практичних проблем рiчкової гiдравлiки [5],
але його використання обмежене випадками, коли
i потiк, i донна поверхня змiнюються дуже повiль-
но.
При моделюваннi нестацiонарних проблем, коли
змiни на доннiй поверхнi можуть впливати на по-
тiк, бiльш перспективним є одночасне розв’язання
гiдродинамiчних i морфологiчного рiвнянь. В ро-
ботi [4] показано, що при зв’язаному iнтегруван-
нi гiдроморфодинамiчної системи чисельнi апро-
ксимацiї є достатньо стiйкими, а результати мо-
делювання - адекватними до фiзичних процесiв.
Проблема полягає в тому, що в бiльшостi випад-
кiв залежностi потоку зважених частинок вiд па-
раметрiв водної течiї мають складний емпiричний
характер, iнколи вони навiть не можуть бути за-
писанi в аналiтичному виглядi. Це значно ускла-
днює знаходження швидкостей поширення збу-
рень у процесах, що розглядаються. Хоча останнiм
часом зв’язаний пiдхiд активно розвивався [6, 7],
одержанi розв’язки стосуються самої простої сте-
пеневої залежностi потужностi потоку наносiв вiд
швидкостi водної течiї.
Зв’язування процесiв може бути частковим, ко-
ли на одному i тому ж кроцi по часовi система
розщеплюється на два оператори, якi виконуються
послiдовно. Цей пiдхiд є бiльш простим, нiж пов-
нiстю зв’язаний, оскiльки характеристичнi швид-
костi у гiдродинамiчнiй i у морфологiчнiй моде-
лях оцiнюються окремо. З iншого боку, вiн здатен
описувати взаємний вплив процесiв, оскiльки ре-
зультати гiдродинамiчного моделювання викори-
стовуються як вхiднi данi при iнтегруваннi рiвня-
ння еволюцiї поверхнi i навпаки. Стiйкiсть таких
алгоритмiв залежить вiд типу i порядку чисельних
схем, якi застосовуються для iнтегрування гiдро-
морфодинамiчних рiвнянь [4, 8, 9].
В цiй роботi розвинений чисельний алгоритм,
який поєднує центральну схему другого поряд-
ку Курганова-Ноелля-Петрової (KNP-scheme) [10,
11] для iнтегрування гiдродинамiчних рiвнянь зi
схемою "проти течiї" п’ятого порядку, за якою
розраховується еволюцiя донної поверхнi. Остан-
ня ґрунтується на зваженiй iстотно неосцiлюючiй
апроксимацiї (WENO-approximation) потоку нано-
сiв, запропонованiй в [12]. Розглянутi одновимiрнi
реалiзацiї цього алгоритму для незв’язаної i час-
тково зв’язаної систем гiдроморфодинамiчних рiв-
нянь. Виконанi тестовi розрахунки еволюцiї пi-
щаного пагорбу, розташованого на днi прямоку-
тного каналу. Порiвняння одержаних результатiв
з вiдомим аналiтичним розв’язком [4] показують,
що обидва пiдходи, незв’язаний i частково зв’яза-
ний, забезпечують прийнятну точнiсть, але остан-
нiй бiльше пiдходить для довгострокових розра-
хункiв, оскiльки майже повнiстю виключає хибнi
коливання шуканих розв’язкiв. Показано, що то-
чнiсть розрахункiв визначається порядком чисель-
ної схеми по простору i мало залежить вiд способу
iнтегрування по часовi, що iстотно спрощує засто-
сування цього алгоритму до моделювання практи-
чних проблем.
Проведенi довгостроковi розрахунки, в яких
отриманi часовi i просторовi характеристики ево-
люцiї глибоких ям (кар’єрiв) та пiщаних пагорбiв
на рiчковому днi. Вони показують, що в процесi
розмиву не лише змiнюється поперечний профiль
нерiвностi, а й вiдбувається її перенос за течiєю.
Зробленi оцiнки є важливими для прогнозуван-
ня впливу крупномасштабних донних утворень на
еволюцiю рiчкових течiй, розмивнi процеси побли-
зу берегiв, стiйкiсть гiдротехнiчних споруд, тощо.
1. ТЕОРЕТИЧНА МОДЕЛЬ
Течiї у природних водоймах, де вертикальний
розмiр є набагато меншим за горизонтальнi, за-
звичай описуються моделлю мiлкої води. Її рiв-
няння виводяться iз загальних законiв динамiки
рiдини iнтегруванням по вертикальнiй координатi
за умови гiдростатичного розподiлу тиску. Якщо
вважати також, що параметри течiї не змiнюються
вздовж поперечного перерiзу русла, можна засто-
совувати одновимiрну модель, яка у випадку пря-
мокутного каналу зводиться до класичної системи
рiвнянь Сен-Венана [13]:
∂h
∂t
+
∂q
∂x
= 0, (1)
∂q
∂t
+
∂
∂x
(
q2
h
+
1
2
gh2
)
= −gh
dzb
dx
− ghSf , (2)
де x – горизонтальна змiнна вздовж вiсi каналу;
t – час; g – прискорення вiльного падiння;
h(x, t) – глибина води; u(x, t) – швидкiсть течiї;
q = hu – витрата води; zb(x, t) – функцiя, яка
задає форму дна; Sf – нахил тертя.
22 I. М. Горбань
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
Рис. 1. Схема мiлкої води в каналi з розмивним дном
Рiвняння (1), (2) є осередненими по вертикаль-
нiй координатi рiвняннями збереження маси i мо-
менту рiдини вiдповiдно. В цiй постановцi вважа-
ється, що конфiгурацiя дна у каналi змiнюється з
часом (рис. 1). Еволюцiя дна описується рiвнян-
ням Екснера, яке випливає iз закону збереження
маси твердої речовини i має наступний вигляд:
(1 − pm)
∂zb
∂t
+
∂qb
∂x
= 0, (3)
де pm – пористiсть донного ґрунту. Функцiя
qb(h, q) представляє собою потужнiсть потоку зва-
жених частинок, якi перемiщуються у придонно-
му шарi. Вона залежить вiд гiдродинамiчних па-
раметрiв течiї та фiзичних властивостей наносiв i
має емпiричний характер.
Рiвняння балансу наносiв (3) описує перенос
твердої речовини лише в межах придонної обла-
стi, тодi як перемiщеннями домiшок у суспензiйно-
му шарi, а також їхнiм випаданням на дно нехту-
ється. Цей пiдхiд є правомiрним, оскiльки в рiв-
нинних водоймах у верхньому шарi переносяться,
головним чином, мул i глина, якi мають слабку
здатнiсть до осiдання. З iншого боку, пiщанi ча-
стинки, якi формують дно, майже не пiднiмаю-
ться вгору через те, що вертикальна складова тур-
булентної швидкостi у такому потоцi є, зазвичай,
меншою за швидкiсть осiдання частинки.
Для замикання системи рiвнянь (1)–(3) вико-
ристовуються емпiричнi формули, якi виражають
тертя Sf i потужнiсть потоку зважених частинок
qb через гiдродинамiчнi змiннi. Величина Sf за-
звичай моделюється законами Дарсi-Вейсбаха або
Маннiнга:
Sf =
fu|u|
8gRh
, (4)
Sf =
n2u|u|
R
4/3
h
, (5)
де f, n – коефiцiєнти Дарсi-Вейсбаха i Маннiнга
вiдповiдно. Їхнi значення для конкретних водойм
вираховуються за допомогою натурних експери-
ментiв.
Величина Rh в формулах (4), (5) є гiдравлi-
чним радiусом, який представляє собою вiдноше-
ння площини поперечного перерiзу потоку до змо-
ченого периметру. У випадку широкого прямоку-
тного каналу, що тут розглядається, Rh ∼ h.
Оцiнки, якi iснують для визначення потужно-
стi потоку наносiв qb, як правило, сформульова-
нi для гранульованих зв’язаних частинок i мають
вигляд iмовiрнiсних або детермiнiстичних законiв.
Покладаючи, що пульсацiї швидкостi в потоцi роз-
подiленi за нормальним законом, Ейнштейн [14]
одержав залежнiсть потоку зважених частинок вiд
iмовiрностi втягнення такої частинки у рух турбу-
лентною течiєю. Детермiнiстичнi моделi зв’язують
транспорт наносiв з фiзичними характеристика-
ми частинок i гiдродинамiчними параметрами те-
чiї. Найбiльш простою є степенева формула Грас-
са [15], в якiй потужнiсть потоку наносiв безпосе-
редньо виражається через гiдродинамiчнi змiннi h
i q:
qb = Ag
q
h
∣∣∣
q
h
∣∣∣
mg−1
, 1 ≤ mg ≤ 4. (6)
Параметр Ag [м2/c] залежить вiд розмiрiв ча-
стинок, що переносяться, i в’язкостi рiдини. Вiн
визначається експериментальним шляхом. Малi
значення Ag (Ag ∼ 0) означають слабку взаємо-
дiю мiж течiєю води i транспортом наносiв. Якщо
Ag ∼ 1, має мiсце сильна взаємодiя.
Зазвичай в формулi (6) покладається mg = 3.
Тодi вона набуває вигляду:
qb = Agu3. (7)
Аналогiчна формула отримана в IГМ НАН Укра-
їни спецiально для умов рiчки Днiпро [2]:
qb = D
u6
h
, (8)
де D – емпiричний коефiцiєнт, D ≈ 2·10−5 [с5/м3].
Бiльш практичнi оцiнки потоку наносiв qb ґрун-
туються на положеннi, що зважена частинка по-
чинає рухатися, коли дотичнi напруження на днi
τb перевищують задане критичне значення τb cr. В
ламiнарному потоцi зв’язок мiж дотичними напру-
женнями τb i силою тертя Sf на доннiй поверхнi
має вигляд:
τb = ρgh|Sf |, (9)
де ρ – густина води; Sf знаходиться за емпiри-
чними законами (4) або (5).
Зазвичай розглядають безрозмiрнi дотичнi на-
пруження τ∗
b i τ∗
b cr:
τ∗
b =
τb
(γs − γ)ds
, τ∗
b cr =
τb cr
(γs − γ)ds
, (10)
I. М. Горбань 23
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
де γ = ρg, γs = ρsg – питома вага води i наносiв
вiдповiдно, ρs – густина наносiв; ds – середнiй
розмiр зважених частинок.
Порогове значення дотичного напруження τ∗
b cr ,
при якому частинка, що лежить на днi, втягується
у рух, залежить вiд параметрiв течiї i властиво-
стей донних вiдкладень. Результати багатьох екс-
периментальних дослiджень по визначенню τ∗
b cr
узагальнив Шильдс, який одержав емпiричну за-
лежнiсть τ∗
b cr вiд числа Рейнольдса придонного
потоку, вiдому, як дiаграма Шiльдса [3]:
τ∗
b cr = f(Re), (11)
де Re = U∗ds/ν , U∗ – дотична швидкiсть потоку
води поблизу дна, ν – кiнематична в’язкiсть рiди-
ни.
У возрахунках використовується безрозмiрна
потужнiсть потоку наносiв:
Φ =
qb√(ρs
ρ
− 1
)
gd3
s
. (12)
Для її визначення розвиненi рiзнi емпiричнi моде-
лi, якi застосовуються в залежностi вiд фiзичних
властивостей наносiв та параметрiв водної течiї.
• Модель Мейера-Петера i Мюллера (1948) [16]:
Φ = 8
(
τ∗
b − τ∗
b cr
)3/2
, (13)
де τ∗
b cr = 0.047.
Ця формула є однiєю з найбiльш вiдомих. Во-
на застосовується, якщо середнiй розмiр зва-
жених частинок d50 = 0.4 мм ÷ 29 мм, а нахил
рiчкового дна не перевищує 2%.
• Рiвняння Ван Рiня (1984) [17]:
Φ =
0.005
C1.7
D
(ds
h
)0.2√
τ∗
b
(
τ∗
b − τ∗
b cr
)2.4
, (14)
де CD – опiр частинки.
• Модель Нiльсена (1992) [18]:
Φ = 12
√
τ∗
b
(
τ∗
b − τ∗
b cr
)
, (15)
де τ∗
b cr = 0.05.
• Модель Рiбберiнка (1998) [3]:
Φ = 11
(
τ∗
b − τ∗
b cr
)1.65
. (16)
Пористiсть донного ґрунту pm представляє со-
бою вiдношення об’єму пустот у речовинi до її за-
гального об’єму:
pm =
Vv
Vv + Vs
, (17)
де Vv – об’єм пустот; Vs – об’єм твердої речови-
ни. Якщо дно складається з однорiдних гранульо-
ваних зв’язаних частинок, для знаходження pm в
роботi [3] запропоноване наступне емпiричне спiв-
вiдношення:
pm =
1 − 0.525
( ds
ds + 4δ
)3
, ds < 1mm,
0.3 + 0.175e
−0.05(ds−d0)
d0 , ds ≥ 1mm,
(18)
де d0 – вiдносний розмiр частинок, зазвичай, по-
кладається d0 = 1 мм; δ – товщина водної плiвки,
що оточує частинку (δ ∼ 0.0004 мм).
2. ЧИСЕЛЬНА МОДЕЛЬ
У цiй роботi чисельна апроксимацiя рiвнянь (1)–
(3) буде виконуватися в рамках напiвзв’язаного
пiдходу, коли система розщеплюється на два опе-
ратори – гiдродинамiчний i морфологiчний, якi ви-
конуються послiдовно в рамках одного часового
кроку. В цьому разi, кожне з рiвнянь системи iн-
тегрується незалежно вiд iншого найбiльш прида-
тним для цього способом, а взаємний вплив про-
цесiв враховується тим, що результати виконання
одного оператору є вхiдними даними для iншого,
i навпаки. Стосовно процесiв, якi тут розглядаю-
ться, це означає, що рiвняння мiлкої води (1)–(2)
розв’язуються при вiдомiй з попереднього кроку
функцiї дна zb, а рiвняння (3) використовується
для оновлення топографiї дна по вже вiдомим зна-
ченням гiдродинамiчних змiнних.
Чисельний алгоритм, розвинений тут для
розв’язання системи рiвнянь (1)–(3), поєднує
центральну схему другого порядку Курганова-
Ноелля-Петрової (KNP-scheme) [10, 11] для iнте-
грування гiдродинамiчних рiвнянь зi схемою "про-
ти течiї" п’ятого порядку, за якою розраховує-
ться еволюцiя донної поверхнi. Остання ґрунтує-
ться на зваженiй iстотно неосцiлюючiй апроксима-
цiї (WENO-approximation) потоку наносiв, запро-
понованiй в [12].
2.1. Чисельна апроксимацiя рiвнянь мiлкої води
Головна складнiсть при побудовi чисельних
апроксимацiй для рiвнянь мiлкої води полягає в
тому, що вони мають гiперболiчну природу, вна-
слiдок чого їхнi розв’язки можуть мiстити роз-
риви. Взаємодiя потоку з сильно нерiвним дном,
яке є предметом дослiдження цiєї роботи, немину-
че призводить до утворення складної конфiгурацiї
вiльної поверхнi, включно зi стрибками i хвилями
розрiдження. Коректна чисельна схема повинна не
24 I. М. Горбань
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
лише описувати такi розв’язки без хибних осциля-
цiй, а й зберiгати їх протягом довгого часу. Одним
iз самих простих i ефективних пiдходiв, якi вiдпо-
вiдають перелiченим вимогам, є центральна схема
"проти течiї" , обґрунтована в роботах [10, 11, 19].
Вона належить до проекцiйно-еволюцiйних мето-
дiв типу Годунова [20], якi мають високу роздiль-
нiсть та можуть знаходити розривнi розв’язки. До-
кладне описання напiвдискретної центральної схе-
ми "проти течiї" представлене у роботi [21], тому
тут наведемо лише загальнi її особливостi.
Система рiвнянь мiлкої води (1), (2) є неоднорi-
дним законом збереження, який у векторному ви-
глядi може бути записаний наступним чином:
Ut + Fx
(
U
)
= S
(
U(x, t), x, t
)
, (19)
де U – вектор консервативних змiнних; F – ве-
ктор потокiв цих змiнних; S – вектор-функцiя,
яка задає джерела втрат консервативних змiнних:
U =
(
h
q
)
, F =
q
q2 +
1
2
gh2
,
S =
0
−gh
(dzb
dx
+ Sf
)
.
Метод скiнченних об’ємiв Годунова для розв’я-
зання рiвняння (19) ґрунтується на його iнте-
гральнiй формi:
dU(x, t)
dt
= −
[
F
(
U(x + ∆x/2, t)
)
−
−F
(
U(x − ∆x/2, t)
)]
+ S
(
U(x, t), x, t
)
,
(20)
де U(x, t) i S
(
U(x, t), x, t
)
– осередненi значення
консервативних змiнних i функцiї джерела втрат
по лiнiйному об’єму I(x) = {ξ : |ξ − x| ≤ ∆x/2}:
U(x, t) ≈ 1
∆x
∫
I(x)
U (ξ, t)dξ,
S
(
U(x, t), x, t
)
≈ 1
∆x
∫
I(x)
S (U(ξ, t), ξ, t)dξ.
(21)
Iнтегральний закон (20) є бiльш загальним, нiж
його диференцiйна форма (19), оскiльки дозволяє
описати розриви у розв’язках, тодi як похiднi в
точках розриву не визначаються.
На однорiднiй розрахунковiй сiтцi в просторi
(x− t) : tn = n∆t, xj = j∆x з рiвняння (20) випли-
ває наступна чисельна схема вiдносно осереднених
консервативних змiнних:
U
n+1
j = U
n
j − ∆t
∆x
[
F̃ n
j+1
2
− F̃ n
j−1
2
]
+ S
n
j , (22)
де функцiя F̃ n
j±1
2
апроксимує осередненi по часовi
потоки змiнних через границi контрольного об’є-
му:
F̃ n
j±1
2
≈
tn+1∫
tn
F
(
U(xj±1
2
, t)
)
dt. (23)
В класичних upwind схемах для розв’язання гi-
перболiчних рiвнянь чисельний потiк F̃ n
j± 1
2
оцiню-
ється за допомогою рiманiвських алгоритмiв (сол-
верiв) [22]. Рiманiвськi солвери, якi обчислюють
значення функцiї вiд розривної змiнної, будуються
на основi характеристичної iнформацiї про вели-
чину i напрямки поширення збурень, що дозволяє
максимально точно описувати розриви у розв’яз-
ках. Центральнi алгоритми ґрунтуються на вiдо-
мiй схемi Лакса-Фрiдрiхса [23]. Вони є значно про-
стiшими та унiверсальнiшими, але протягом дов-
гого часу їхнiм iстотним недолiком була велика
чисельна дифузiя. В роботi [19] розроблений но-
вий клас центральних схем, де iнтегрування за
часом пропонується виконувати окремо по тiй ча-
стинi контрольного об’єму, де змiнна неперервна,
i по областi розриву, ширина якої оцiнюється за
допомогою локальних швидкостей поширення збу-
рень. Застосування характеристичної iнформацiї
дозволяє iстотно зменшити чисельну дифузiю цен-
тральних алгоритмiв, наблизивши їх до схем "про-
ти течiї".
В напiвдискретних центральних схемах iнтегру-
вання по простору i по часовi виконується послi-
довно, отже, вони зводяться до системи звичай-
них диференцiальних рiвнянь вiдносно осередне-
них консервативних змiнних:
dUj(t)
dt
= −
Hj+1
2
(t) − Hj−1
2
(t)
∆x
+ Sj(t), (24)
де Hj±1
2
(t) – чисельнi потоки змiнних крiзь грани-
цi j-го контрольного об’єму.
Вигляд функцiї Hj+ 1
2
(t) залежить вiд способу
апроксимацiї (реконструкцiї) змiнних на грани-
цях елементарних об’ємiв. Схема другого поряд-
ку Курганова-Ноелля-Петрової [10] ґрунтується
на лiнiйнiй реконструкцiї змiнних:
Ũ(x, t) = Uj(t) + (Ux)j(x − xj), x ∈ Ij . (25)
Для мiнiмiзацiї коливань реконструйованих
змiнних при визначеннi похiдних у (25) використо-
вуються спецiальнi нелiнiйнi функцiї вiд локаль-
них градiєнтiв, якi називаються лiмiтерами. В цiй
I. М. Горбань 25
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
роботi застосовується однопараметрична minmod-
функцiя [24]:
(Ux)j = minmod
(
θ
U j+1 − U j
∆x
,
Uj+1 − U j−1
2∆x
, θ
U j − U j−1
∆x
)
,
(26)
де θ ∈ [1, 2]. Зауважимо, що коливання змiнних
будуть меншими при малих значеннях параметру
θ iз зазначеного дiапазону, але дисипацiя чисель-
ної схеми при цьому збiльшується.
Застосування лiмiтеру робить цю чисельну схе-
му нелiнiйною, отже, її розв’язки мають принайм-
нi другий порядок точностi в областях монотон-
ностi змiнних, крiм того, усуваються коливання
шуканих величин поблизу великих градiєнтiв або
розривiв. Завдяки цьому, схема, що розглядається,
належить до класу TVD-методiв (Total Variation
Diminishing), для яких помилка дискретизацiї з
часом не збiльшується [25].
Побудова чисельного потоку для схеми (24)
докладно описана у роботi [10]. Тут наведемо
остаточний вигляд функцiї H(t):
Hj+1
2
(t) =
=
a+
j+ 1
2
F (U−
j+ 1
2
) − a−
j+ 1
2
F (U+
j+ 1
2
)
a+
j+ 1
2
− a−
j+ 1
2
+
+a+
j+ 1
2
a−
j+ 1
2
(
U+
j+ 1
2
− U−
j+ 1
2
a+
j+ 1
2
− a−
j+ 1
2
− qj+ 1
2
)
,
(27)
де U±
j+ 1
2
, a±
j+1
2
– значення консервативних змiн-
них та екстремальних хвильових швидкостей злiва
i справа вiд границi xj+ 1
2
вiдповiдно.
Швидкостi a±
j+ 1
2
описують ширину областi роз-
риву змiнної U±
j+ 1
2
на границi елементарного об’є-
му. Вони знаходяться наступним чином:
a−
j+ 1
2
=
= min
{
λ1
(
J(U−
j+ 1
2
)
)
, λ1
(
J(U+
j+ 1
2
)
)
, 0
}
,
a+
j+ 1
2
=
= max
{
λ2
(
J(U−
j+ 1
2
)
)
, λ2
(
J(U+
j+ 1
2
)
)
, 0
}
,
(28)
де λ1, λ2 – власнi значення Якобiану J = ∂F/∂U
рiвняння (19), причому λ1 = u−
√
gh, λ2 = u+
√
gh.
Величина qj+1
2
в (27) вводиться для зниження си-
стемної дифузiї на проекцiйному кроцi. Її вигляд
можна знайти в [11].
Iнтегрування по часовi в проекцiйно-
еволюцiйних методах, до яких належить чи-
сельний алгоритм, що розглядається, зазвичай
виконується за явною схемою. Її стiйкiсть забезпе-
чується використанням багатокрокових SSP-схем
високих порядкiв (Strong Stability Preserving) [26],
якi будуються на основi методу Рунге-Кутта.
Крок по часовi в цих схемах повинен вiдповiдати
критерiю Куранта-Фрiдрiхса-Левi:
a∆t
∆x
≤ CFL (29)
при CFL = 0.5. Величина a в (29) є максималь-
ною швидкiстю поширення збурень у розрахунко-
вiй областi на заданому часовому кроцi:
a = maxj
{
max{a+
j+ 1
2
, a−
j+1
2
}
}
.
Важливою вимогою до чисельних алгоритмiв,
якi створюються для розв’язання неоднорiдних гi-
перболiчних рiвнянь, подiбних до (19), є збереже-
ння стацiонарних розв’язкiв з порядком точностi
не нижче другого. Тобто, при моделюваннi ста-
цiонарних проблем повинен виконуватися баланс
мiж лiвою i правою частинами рiвнянь. В схемi,
що розглядається, це забезпечується застосуван-
ням кусково-лiнiйної апроксимацiї дна i спецiаль-
них квадратур для дискретизацiї правої частини
рiвняння [27].
Iнша складнiсть полягає в тому, щоб забезпечи-
ти коректне моделювання гiдродинамiчних змiн-
них, коли глибина води наближається до нуля,
зокрема, при поширеннi фронтiв змочення поверх-
нi. Тодi коливання змiнної поблизу розривiв мо-
жуть зупинити розрахунок через те, що вона може
стати негативною. В цiй схемi додатнiсть глибини
потоку у всiй областi течiї досягається за допомо-
гою спецiальної процедури корекцiї змiнних, вве-
деної в [27].
Виконанi в [21] докладнi тестування центральної
схеми проти течiї на класичних гiдравлiчних при-
кладах i даних лабораторних вимiрювань в одно-
i двовимiрном випадках свiдчать про те, що пред-
ставлена чисельна модель:
– добре описує малi збурення вiльної поверхнi,
– дозволяє прогнозувати розповсюдження фронтiв
змочення на днi (wet/dry fronts),
– здатна зберiгати стацiонарнi гiдравлiчнi розв’яз-
ки,
– адекватно описує генерацiю хвиль при руйнуван-
нi дамб i гребель,
– дозволяє описувати гiдравлiчнi течiї в каналах зi
складною топографiєю дна, у тому числi там, де є
дiлянки незмоченної донної поверхнi.
26 I. М. Горбань
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
2.2. Схема високого порядку для рiвняння
балансу наносiв
Рiвняння еволюцiї донної поверхнi (3) часто за-
писують у наступному виглядi:
∂zb
∂t
+ ξ
∂qb
∂x
= 0, (30)
де ξ = 1/(1 − pm).
Оскiльки поле течiї, а за ним i потужнiсть пото-
ку наносiв залежать вiд форми донної поверхнi,
рiвняння (30) представляє собою нелiнiйний закон
збереження вiдносно змiнної zb. Для його розв’я-
зання повиннi застосовуватися методи, якi мо-
жуть знаходити розривнi розв’язки (shock capturi-
ng schemes). Традицiйним алгоритмам, якi ви-
користовують центральнi схеми Лакса-Вендорфа
або Мак-Кормака, властива наявнiсть високоча-
стотних хибних осциляцiй у розв’язках. В роботi
[4] показано, що цi коливання є наслiдком нелiнiй-
ної природи морфологiчного рiвняння (30). Для
згладжування розв’язкiв у цих алгоритмах засто-
совуються рiзнi реґуляризацiйнi процедури, на-
приклад, введення штучної в’язкостi [28] або фiль-
трування [29]. Схеми проти течiї, якi будуються на
основi рiманiвських солверiв [4], є бiльш стiйкими,
нiж центральнi, але мають низький порядок то-
чностi, тодi як акуратне моделювання морфологi-
чних процесiв потребує схем високого порядку.
В роботi [12] представлений новий алгоритм для
розв’язання еволюцiйного рiвняння (30), в яко-
му висока точнiсть результатiв досягається вико-
ристанням зваженої iстотно неосцилюючої апро-
ксимацiї потоку наносiв (WENO-approximation).
WENO-схеми вперше були введенi в 1994 роцi в ро-
ботi [30] i з тих пiр активно залучаються до чисель-
ного iнтегрування диференцiальних рiвнянь у ча-
стинних похiдних, у тому числi гiперболiчного ти-
пу. Їхньою перевагою є високий порядок точностi
в областях монотонностi шуканих величин i стiй-
ка, неосцилююча поведiнка одержуваних розв’яз-
кiв поблизу розривiв. Проведенi в [12] розрахунки
тестової морфологiчної задачi про розмив пiща-
ного пагорбу виявили iстотнi переваги запропоно-
ваної WENO-схеми над класичними алгоритмами
щодо моделювання деформацiй донної поверхнi.
В цiй роботi WENO-схема з [12] буде iнтегрована
в загальну чисельну модель, яка виконує зв’язанi
розрахунки гiдродинамiчних i морфологiчних па-
раметрiв потоку.
Сучаснi shock-capturing алгоритми для розв’я-
зання рiвняння (30) ґрунтуються на його хвильо-
вiй апроксимацiї i використовують фазовi швид-
костi поширення збурень на доннiй поверхнi [28].
Як вже вiдмiчалося, потужнiсть потоку наносiв qb
залежить вiд багатьох факторiв - швидкостi течiї,
глибини води, характеристик хвиль, фiзичних вла-
стивостей частинок, тощо. Але, якщо розглядати
морфологiчний процес на окремому кроцi по ча-
совi, коли рiвень води, гiдродинамiчнi параметри
течiї i властивостi наносiв фiксованi, то qb можна
вважати функцiєю лише вiд рiвня дна zb, i тодi
виконується:
∂qb
∂x
=
∂qb
∂zb
∂zb
∂x
. (31)
З урахуванням (31) рiвняння (30) приймає вигляд:
∂zb
∂t
+ C(zb)
∂zb
∂x
= 0, (32)
де C(zb) = ξ
∂qb
∂zb
– фазова швидкiсть поширення
збурень на доннiй поверхнi.
Очевидно, що (32) є рiвнянням адвекцiї вiдно-
сно змiнної zb. З нього випливає, що морфологi-
чна еволюцiя вiдбувається як нелiнiйне пошире-
ння деформацiй дна. Перенос форми дна в нега-
тивному i позитивному напрямках характеризує-
ться фазовими швидкостями C− = min{C, 0} i
C+ = max{C, 0} вiдповiдно, так що
C(z) = C−(z) + C+(z). (33)
Транспорт наносiв також складається з двох ча-
стин:
qb = qb− + qb+, (34)
де
q−b = (1 − pm)
zb∫
0
C−(z)dz,
q+
b = (1 − pm)
zb∫
0
C+(z)dz.
(35)
Подiбно до методу скiнченних об’ємiв, дискре-
тизацiя рiвняння (30) у схемi, що розглядається,
виконується на сiтцi, яка складається з вiдрiз-
кiв [xi−1
2
, xi+1
2
], так що точки xi = i∆x, в яких
оновлюється дно, лежать посерединi розрахунко-
вих елементiв. Для обчислення ∂qb/∂x в цих то-
чках застосовується схема першого порядку:
(
∂qb
∂x
)
i
=
q̂bi+1
2
− q̂bi−1
2
∆x
. (36)
Величини q̂bi±1
2
в (36) є чисельними апроксимацi-
ями потужностi потоку наносiв qb на кiнцях еле-
ментарного об’єму. Головна задача при побудовi
чисельної схеми полягає в тому, щоб вона як мо-
жна краще оцiнювала q̂bi±1
2
, оскiльки вiд цього за-
лежить якiсть одержуваних результатiв.
I. М. Горбань 27
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
В роботi [12] пропонується розщеплювати чи-
сельнi потоки на границях контрольних елемен-
тiв на лiвостороннiй i правостороннiй, якi потiм
обчислюються з урахуванням напрямку пошире-
ння збурень на доннiй поверхнi. Завдяки цьому,
одержана схема має природу "проти течiї" i, як
всi аналогiчнi пiдходи, здатна забезпечити високу
стiйкiсть розрахункiв. Вiдповiдно до (34) в точцi
xi+1
2
маємо:
q̂bi+1
2
= q̂−
bi+1
2
+ q̂+
bi+1
2
. (37)
Одностороннi потоки q̂±
bi+ 1
2
обчислюються, вихо-
дячи зi значень qb в сусiднiх до xi+ 1
2
точках. Чим
бiльше точок застосовується для оцiнки, тим ви-
щим є порядок чисельної схеми.
Зважена iстотно неосцiлююча upwind апрокси-
мацiя для лiвостороннього потоку наносiв q̂−
bi+1
2задається виразами:
q̂−
bi+1
2
=
ω1q
(1)
i+ 1
2
+ ω2q
(2)
i+1
2
+ ω3q
(3)
i+1
2
, Ci+1
2
≥ 0,
0, Ci+ 1
2
< 0,
(38)
де
q
(1)
i+ 1
2
=
1
3
qbi−2 −
7
6
qbi−1 +
11
6
qbi,
q
(2)
i+ 1
2
= −1
6
qbi−1 +
5
6
qbi +
1
3
qbi+1,
q
(3)
i+ 1
2
=
1
3
qbi +
5
6
qbi+1 −
1
6
qbi+2.
(39)
Ця апроксимацiя використовує три точки, якi ле-
жать злiва вiд границi xi+ 1
2
(з iндексами i−2, i−1,
i), отже, має 3-й порядок точностi. Подальше пiд-
вищення порядку схеми досягається за рахунок ва-
гових коефiцiєнтiв. В роботi [31] описаний спосiб
розрахунку коефiцiєнтiв ω1, ω2, ω3, при яких схема
(38) має 5-й порядок. Наведемо тут кiнцевi резуль-
тати цього дослiдження:
ω1 =
α1
Sα
, ω2 =
α2
Sα
, ω3 =
α3
Sα
, (40)
де Sα = α1 + α2 + α3,
α1 =
0.1
(S1 + ε)2
, α2 =
0.6
(S2 + ε)2
,
α3 =
0.3
(S3 + ε)2
.
(41)
Величина ε в формулах (41) є малим числом, яке
вводиться, щоб виключити дiлення на нуль. Дис-
кретнi функцiї S1, S2, S3 називаються мiрою згла-
джування i виражаються через потiк наносiв у су-
сiднiх до границi xi+ 1
2
точках. Їхнiй вигляд можна
знайти в [12].
Правостороннiй потiк q̂+
bi+1
2
розраховується за
аналогiчною технологiєю:
q̂+
bi+1
2
=
ω4q
(4)
i+1
2
+ ω5q
(5)
i+1
2
+ ω6q
(6)
i+ 1
2
, Ci+1
2
< 0,
0, Ci+ 1
2
≥ 0,
(42)
q
(4)
i+1
2
= −1
6
qbi−1 +
5
6
qbi +
1
3
qbi+1,
q
(5)
i+1
2
=
1
3
qbi +
5
6
qbi+1 −
1
6
qbi+2,
q
(6)
i+1
2
=
11
6
qbi+1 −
7
6
qbi+2 +
1
3
qbi+3.
(43)
Ваговi коефiцiєнти ω4, ω5, ω6 обчислюються за
формулами (40), де
α1 =
0.3
(S4 + ε)2
, α2 =
0.6
(S5 + ε)2
,
α3 =
0.1
(S6 + ε)2
.
(44)
Вигляд функцiй S4, S5, S6 є в роботi [12].
Аналогiчно оцiнюється потiк на границi xi− 1
2
,
лише всi iндекси у формулах зсуваються на оди-
ницю влiво.
Для обчислення хвильових швидкостей на
границях елементiв сiтки в класичних алгори-
тмах ”проти течiї” використовується центральна
кiнцево-рiзницева формула [4]:
Ci+ 1
2
= ξ
qbi+1 − qbi
zbi+1 − zbi
, (45)
яка втрачає чиннiсть, коли zbi+1 ≈ zbi. Перева-
гою представленої в цiй роботi схеми є те, що вона
потребує iнформацiї лише про знак швидкостi, i
тодi замiсть дiлення в (45) може застосовуватися
множення (qbi+1 − qbi)(zbi+1 − zbi), яке виключає
невизначенiсть при (zbi+1 − zbi) ≈ 0.
Якщо для дискретизацiї рiвняння (30) по часовi
застосувати метод Ейлера, то остаточна чисель-
на схема для розрахунку еволюцiї донної поверхнi
набуває вигляду:
zn+1
bi = zn
bi −
∆tm
∆x
ξ
(
q̂bi+1
2
− q̂bi−1
2
)
, (46)
де ∆tm – крок дискретизацiї за часом морфологi-
чного процесу.
В роботi [12] показано, що застосування бага-
токрокових схем типу Рунге-Кутта до моделюва-
ння динамiки дна слабо впливає на результати, в
той самий час, воно потребує перерахунку гiдроди-
намiчних параметрiв течiї на промiжних кроках.
Такий пiдхiд може виявитися занадто затратним
28 I. М. Горбань
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
при виконаннi практичних розрахункiв. Представ-
лений тут алгоритм, який поєднує WENO-схему з
методом Ейлера, є оптимальним щодо точностi ре-
зультатiв i простоти реалiзацiї. Зазначимо також,
оскiльки схема (46) є явною, вона потребує вико-
нання стандартного критерiю стiйкостi Куранта-
Фрiдрiхса-Левi:
c
∆tm
∆x
≤ 1, (47)
де c = maxj
{
|Ci+1
2
|
}
.
3. ТЕСТОВI РОЗРАХУНКИ
В цьому роздiлi оцiнюється здатнiсть представ-
леної чисельної моделi розраховувати гiдродинамi-
чнi та морфологiчнi змiни у вiдкритих водоймах.
Для того щоб виявити, який з двох пiдходiв, ста-
цiонарний чи нестацiонарний, забезпечує вищу то-
чнiсть результатiв, гiдроморфодинамiчна пробле-
ма (1)–(3) розглядається у двох постановках.
У першiй постановцi, позначимо її А1, вважа-
ється, що морфологiчнi змiни на доннiй поверхнi
вiдбуваються так повiльно, що не впливають на
течiю. Тодi рiвняння мiлкої води (1)–(2) розв’язу-
ються при незмiнному днi до тих пiр, поки течiя не
досягне стану рiвноваги, що вiдповiдає виконанню
наступного критерiю:∣∣wn+1
i − wn
i
∣∣ ≤ tol, (48)
де wi – рiвень вiльної поверхнi у вузлах розрахун-
кової сiтки; tol – задане мале значення. Пiсля цьо-
го iнтегрується еволюцiйне рiвняння (3) для дна
за умови, що всi параметри водної течiї фiксова-
нi. Очевидно, що в цьому формулюваннi хвильовi
швидкостi поширення збурень у гiдродинамiчному
i морфологiчному процесах iстотно вiдрiзняються,
отже, морфологiчний крок по часовi ∆tm є наба-
гато бiльшим, нiж гiдродинамiчний ∆t. Загальний
час процесу визначається за морфологiчним кро-
ком.
В нестацiонарнiй постановцi, позначеної як А2,
нiяких умов на течiю i еволюцiю дна не наклада-
ється, а гiдродинамiчнi i морфологiчне рiвняння
розв’язуються в межах одного й того самого кро-
ку по часовi ∆t.
Тестовi розрахунки, виконанi в цiй роботi, вклю-
чають моделювання вiдомих морфологiчних про-
блем, якi мають аналiтичнi розв’язки [4, 12]. Крiм
того, що одержанi чисельнi результати порiвню-
ються з точними розв’язками, у всiх морфологi-
чних розрахунках перевiряється збереження об’є-
му твердої речовини. Її загальний об’єм у заданiй
областi [B1, B2] вираховується за формулою:
V (t) =
B2∫
B1
(zb(x, t) − zb∞) dx, (49)
де zb∞ – умовний рiвень незбуреного дна. Вiдносна
помилка розрахункiв обчислюється як
Error =
V (t) − V (0)
V (0)
. (50)
3.1. Розмив пiщаного пагорбу у стацiонарному
потоцi
Цей класичний приклад застосовується для
апробацiї чисельних схем, якi моделюють транс-
порт наносiв у придоннiй областi [12]. Розрахову-
ється розмив симетричного iзольованого пагорбу,
розташованого на днi каналу, в умовах, коли гi-
дродинамiчнi параметри течiї не змiнюються.
Початкову конфiгурацiю дна задає функцiя:
zb(x, 0) = −h0 + 2 exp [−β(x − xc)] , (51)
де h0 – глибина води над горизонтальним дном;
xc – центр пагорбу; β – параметр, який визначає
крутизну пагорбу. Покладається:
h0 = 6 м, xc = 150 м, β = 0.01.
Довжина розрахункової областi дорiвнює 300 м.
Рiвень та витрата води протягом всього розрахун-
ку вважаються незмiнними:
w(x, t) = 0, q(x, t) = 10 м2/с.
Потужнiсть потоку зважених частинок у розра-
хунках оцiнюється за степеневою формулою Грас-
са (7) при Ag = 0.001 с2/м, пористiсть ґрунту
pm = 0.4. Параметри розрахункової схеми:
∆x = 1 м, ∆t = 0.1 с.
Результати розрахункiв та їхнє порiвняння з то-
чним розв’язком [12] представленi на рис. 2. Тут
можна бачити, що трансформацiї пагорбу внаслi-
док розмиву призводять до формування гострого
фронту седiментних частинок, який рухається за
течiєю. Важливим є те, що чисельна схема описує
змiни конфiгурацiї дна без хибних осциляцiй i не
втрачає стiйкостi протягом достатньо довгого роз-
рахункового промiжку. Одержанi чисельнi резуль-
тати добре узгоджуються з точним розв’язком.
3.2. Розмив пiщаного пагорбу
у нестацiонарному потоцi
На цьому прикладi з роботи [4] тестується чи-
сельнiй алгоритм для розв’язання повної системи
гiдроморфодинамiчних рiвнянь. Тут моделюються
як еволюцiя донної поверхнi, так i зумовленi нею
змiни гiдродинамiчних параметрiв течiї. Розгляда-
ється канал довжиною 1000 м, на днi якого, при
I. М. Горбань 29
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
Рис. 2. Змiнювання розмиву пiщаного пагорбу у
стацiонарному потоцi. Порiвняння чисельних
результатiв з точним розв’язком:
1 – t = 0; 2 – t = 600 c; 3 – t = 2000 c; 4 – t = 6000 c.
300 м ≤ x ≤ 500 м, розташований симетричний
iзольований пагорб висотою 1 м, форма якого опи-
сується функцiєю:
zb(x, 0) = sin2
(π(x − 300)
200
)
.
На початку розрахункiв поверхня води незбуре-
на, а витрата води скрiзь в областi постiйна. По-
кладається:
w(x, 0) = 10 м, q(x, 0) = 10 м2/с.
Просторова сiтка має крок ∆x = 5 м, крок по ча-
совi ∆t в гiдродинамiчних розрахунках обчислює-
ться з критерiю Куранта-Фрiдрiхса-Левi (29) при
CFL = 0.5, ∆tm = 300 с. Зазначимо також, що
на обох границях розрахункової областi накладає-
ться умова вiдкритого каналу, яка означає вiльне
перетiкання води i наносiв через границю.
Потужнiсть седiментного потоку оцiнюється за
степеневою формулою Грасса (7) при Ag = 1 с2/м
та Ag = 0.001 с2/м, що описує сильну i слабку вза-
ємодiю транспорту наносiв з водною течiєю, вiдпо-
вiдно, пористiсть ґрунту pm = 0.4.
Повна гiдроморфодинамiчна проблема розгля-
дається в обох постановках – А1 i А2. Щоб
уникнути рiзкого старту морфологiчного процесу,
спочатку знаходяться параметри стацiонарної те-
чiї, якi вона має при фiксованiй формi дна. Крите-
рiєм досягнення стацiонарного стану є умова (48)
при tol = 10−8. На рис. 3 наведенi початкова кон-
фiгурацiя дна, а також форма вiльної поверхнi
у стацiонарному потоцi над нерозмивним пагор-
бом. З цих результатiв випливає, що над випуклим
дном поверхня води прогинається. Оскiльки па-
горб, що тут розглядається, має незначну висоту у
порiвняннi з глибиною води, то викривлення вiль-
ної поверхнi є дуже слабким, i його не видно на
основному графiку. Залежнiсть, представлена на
Рис. 3. Початкова конфiгурацiя дна i рiвень води в
стацiонарному потоцi над нерозмивним пагорбом
допомiжному рисунку, свiдчить про те, що поверх-
ня води над таким пагорбом прогинається трохи
бiльше, нiж на 1 мм. Здатнiсть застосованої тут
чисельної схеми описувати такi незначнi змiни па-
раметрiв течiї без хибних коливань пiдтверджує її
високу точнiсть щодо розв’язання рiвнянь мiлкої
води.
Рис. 4. Форма розмивного пагорбу при t = 238079 с,
Ag = 0.001 с2/м. Порiвняння чисельних результатiв з
точним розв’язком
Отриманi параметри стацiонарної течiї над не-
розмивним пагорбом покладаються у якостi по-
чаткових умов повної гiдроморфологiчної пробле-
ми, яка розраховує динамiку розмиву i зумовле-
нi нею змiни параметрiв водної течiї. Рис. 4
iлюструє конфiгурацiю розмивного дна zb(x) при
t = 238079 с (≈ 65 годин), яка одержана при слаб-
кiй взаємодiї транспорту наносiв з течiєю води.
Тут порiвнюються результати розрахункiв незв’я-
30 I. М. Горбань
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
заної А1 i зв’язаної А2 задач мiж собою, а та-
кож з наближеним аналiтичним розв’язком, одер-
жаним в роботi [4]. Можна бачити, що чисельна
схема забезпечує акуратне моделювання придон-
ного транспорту наносiв в обох випадках. Це про-
являється у вiдсутностi хибних коливань у чисель-
них результатах та їхньої вiдповiдностi до аналiти-
чного розв’язку.
Разом iз топографiєю дна змiнюються також гi-
дродинамiчнi параметри течiї. Рис. 5,а iлюструє
форму вiльної поверхнi, яка вiдповiдає дну, пред-
ставленому на рис. 4. Можна бачити, що в процесi
розмиву пагорбу генерується перепад глибин (ска-
чок), який рухається разом з донним фронтом се-
дiментних частинок. Викривлений профiль швид-
костi над розмитим пагорбом (рис. 5,b) демонструє
наявнiсть великих градiєнтiв швидкостi над фрон-
тальною частиною нерiвностi. На рис. 4, 5 видно,
що в рамках розглянутого часового перiоду (до
70 годин) результати моделювання за схемами А1
i А2 практично не вiдрiзняються.
Щоб прояснити, чи не втрачається стiйкiсть чи-
сельної схеми з часом, ця тестова задача розра-
ховувалася з тими самими параметрами до t =
150 годин (5400000 с). На рис. 6 показанi конфiгу-
рацiї донної поверхнi, одержанi в незв’язанiй А1
(крива 2) i зв’язанiй А2 (крива 3) постановках.
З них випливає, що оцiнки напрямку i швидкостi
транспорту наносiв в обох пiдходах спiвпадають,
але наявнiсть осциляцiй на кривiй 2 вказує на те,
що застосування стацiонарної моделi з часом може
призвести до втрати стiйкостi чисельної схеми.
Перевага зв’язаної постановки А2 у довгостро-
кових розрахунках пiдтверджується оцiнкою вiд-
носної помилки (50) щодо збереження об’єму твер-
дої речовини в чисельному моделюваннi. Зобра-
женi на рис. 7 залежностi помилки вiд часу по-
казують, що стацiонарний пiдхiд А1 є бiльш ко-
ректним у короткострокових розрахунках, але на
великому перiодi його помилка збiльшується до-
сить швидко, що, очевидно, пов’язане з генерацiєю
хибних коливань у результатах. Нестацiонарна по-
становка А2 забезпечує високу точнiсть результа-
тiв, але на її реалiзацiю потрiбно набагато бiльше
комп’ютерного часу, оскiльки морфологiчна ево-
люцiя розраховується тут з тим самим кроком ∆t,
що й гiдродинамiчна. Зазначимо, що в обох випад-
ках помилка не перевищує порядок 10−4, це свiд-
чить про коректнiсть представленої тут чисельної
схеми щодо iнтегрування законiв збереження.
На рис. 8 показанi динамiка пiщаного пагорбу
i вiдповiднi до неї змiни на поверхнi води, розра-
хованi при сильнiй взаємодiї потоку води i транс-
порту наносiв (Ag = 1 с2/м). Для розв’язання
a
b
Рис. 5. Параметри гiдродинамiчної течiї над
розмивним пагорбом при t = 238079 с:
a – форма вiльної поверхнi; b – швидкiсть течiї
цiєї задачi використовується нестацiонарна поста-
новка А2, яка враховує взаємний вплив процесiв.
Представленi результати демонструють еволюцiю
дна протягом 15 хвилин. Тут можна бачити швид-
ке перемiщення пагорбу за течiєю i формування
фронту седiментних частинок, що супроводжує-
ться утворенням хвилi на поверхнi води. З часом,
амплiтуда цiєї хвилi падає пропорцiйно до зменше-
ння висоти пагорбу. Стiйкiсть розрахункiв i вiдсу-
тнiсть малих хибних коливань у результатах пока-
зують, що застосований чисельний алгоритм зда-
тен моделювати гiдравлiчнi процеси з рiзкою змi-
ною параметрiв.
I. М. Горбань 31
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
Рис. 6. Форма розмивного пагорбу при t = 150 годин:
1 – початкова конфiгурацiя пагорбу; 2, 3 – пагорб,
розрахований в моделях А1 i А2, вiдповiдно
Рис. 7. Вiдносна помилка морфологiчних
розрахункiв:
1 – у стацiонарнiй моделi А1;
2 – у нестацiонарнiй моделi А2
4. МОДЕЛЮВАННЯ ДИНАМIКИ
КРУПНОМАСШТАБНИХ
НЕРIВНОСТЕЙ РIЧКОВОГО ДНА
Деформацiї великих просторових масштабiв на
доннiй поверхнi можуть бути як частиною приро-
дного рельєфу, так i наслiдком iнтенсивної проми-
слової дiяльностi у береговiй зонi. Такi утворен-
ня викликають порушення балансу наносiв i по-
ширення ерозiї в рiчкових руслах, що може при-
звести до пошкоджень та руйнувань iнженерних
систем, якi експлуатуються пiд водою (мостових
опор, трубопроводiв, кабелiв та iн.). Тому дослiд-
ження еволюцiї нерiвностей на рiчковому днi є ва-
жливим для попередження небезпечних екологi-
чних та техногенних явищ в рiчкових басейнах.
Рис. 8. Динамiка розмивного пагорбу та еволюцiя
вiльної поверхнi при сильнiй взаємодiї потоку води
i транспорту наносiв:
1 – t = 200 с; 2 – t = 500 с; 3 – t = 800 с
4.1. Розмив пiщаного пагорбу
Тут дослiджуються механiзм розмиву пiдводно-
го пагорбу, а також зворотнiй вплив цього про-
цесу на гiдродинамiчну течiю. Розглядається ка-
нал довжиною 1000 м, на днi якого розташований
симетричний пiщаний пагорб, локалiзований мiж
300 i 500 м. Головним чинником, який визначає де-
формацiї пагорбу в процесi розмиву, є швидкiсть
потоку води навколо нього. Її величина залежить
не лише вiд витрати води в каналi, а й вiд того,
наскiльки нерiвнiсть перекриває потiк (вiд коефi-
цiєнту запирання). Для дослiдження впливу цiєї
характеристики на розмив будемо розглядати пiд-
воднi пагорби рiзної висоти d = 1 м, d = 1.5 м i
d = 2 м при одних i тих же початкових витратi q i
рiвнi води w в каналi. Покладається:
w(x, 0) = 5 м, q(x, 0) = 10 м2/с.
Потужнiсть потоку седиментних частинок оцi-
нюється за степеневим законом (7) при Ag =
0.005 с2/м, пористiсть ґрунту pm = 0.4. Iнтегрува-
ння рiвнянь (1)–(3) виконується в рамках зв’язаної
постановки А2 при ∆x = 2.5 м, CFL = 0.5.
Спочатку iнтегруються рiвняння мiлкої води до
тих пiр, поки течiя над нерозмивним дном не на-
буде стану рiвноваги. Одержанi розподiли швид-
костi води u вздовж горизонтальної осi каналу x
над пагорбами рiзної висоти d у стацiонарному по-
тоцi представленi на рис. 9. Вони показують, що
залежнiсть швидкостi води вiд висоти пагорбу має
суто нелiнiйний характер, що зумовлене природою
закону збереження (1)–(2). Виконанi за результа-
тами рис. 9 оцiнки числа Фруда (Fr = u/
√
gh)
показують, що над пагорбом висотою 1 м i 1.5 м
32 I. М. Горбань
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
Рис. 9. Залежнiостi швидкостi потоку води u вiд
горизонтальної координати x при рiзнiй висотi
пагорбу:
1 – d = 1 м; 2 – d = 1.5 м; 3 – d = 2 м
ця характеристика змiнюється вiд 0.3 до 0.4, тоб-
то потiк залишається докритичним у всiй обла-
стi. В той самий час, рiзке збiльшення швидкостi
води над пагорбом висотою 2 м супроводжується
переходом течiї у закритичний режим, для якого
Fr ≥ 1.
На рис. 10 суцiльними лiнiями показанi поча-
тковi конфiгурацiї дна i вiдповiднi до них фор-
ми вiльної поверхнi для стацiонарної течiї. Можна
бачити, що у докритичному потоцi поверхня во-
ди над пагорбом плавно прогинається (рис. 10,a).
Амплiтуда прогину складає приблизно 0.5 м. Пе-
рехiд течiї у закритичний режим супроводжується
рiзким перепадом рiвня води над пагорбом (гi-
дравлiчним стрибком), величина якого перевищує
1.5 м (рис. 10,b). В цiй областi локалiзується зна-
чна енергiя води, що може викликати швидкi не-
рiвномiрнi змiни як у потоцi, так i на днi каналу.
Моделювання розмиву пагорбу показує, що пер-
шими у потiк наносiв втягуються частинки пiску,
розташованi на пiдвiтряному боцi пагорбу. Пiсок,
який знаходиться нагорi, переноситься швидше,
нiж той, що прилягає до пiднiжжя нерiвностi. Як
результат, верхня частина пагорбу зсувається на-
перед, i, таким чином, генерується рiзкий фронт
седiментних частинок, який рухається за течiєю
води. У закритичному потоцi (рис. 10,b) розмив
починається iз сильного гiдравлiчного удару об
поверхню пагорбу, внаслiдок чого вiн деформує-
ться. Але надалi поверхня пагорбу згладжується,
i процес протiкає за тим же сценарiєм, що i на
рис. 10,a. По мiрi розмиву пагорб розтягується i
набуває форми довгої нерiвностi з гострим пере-
днiм краєм. Розтягування i вiдповiдне зменшення
висоти пагорбу продовжуються до тих пiр, поки
профiль донної поверхнi не вирiвняється. З рис. 10
a
b
Рис. 10. Динамiка розмиву пiщаного пагорбу i
еволюцiя вiльної поверхнi при рiзнiй висотi пагорбу:
a – d = 1.5 м; b – d = 2 м
випливає, що при однакових гiдродинамiчних умо-
вах на входi у канал нерiвностi бiльшого масштабу
трансформуються швидше, нiж меншi, що, очеви-
дно, зумовлене бiльшою швидкiстю води навколо
великих нерiвностей.
З рис. 10 випливає також, що форма вiльної по-
верхнi в процесi розмиву пагорбу змiнюється вiд-
повiдно до еволюцiї дна. Одночасно з формуван-
ням фронту седиментних частинок на поверхнi во-
ди утворюється хвиля, яка рухається паралельно
перемiщенням твердої речовини на днi каналу.
4.2. Обтiкання i розмив донних заглиблень
В цьому роздiлi представленi результати чи-
сельного моделювання транспорту наносiв навко-
ло донних заглиблень великих просторових мас-
штабiв. Такi заглиблення часто утворюються при
розробцi руслових кар’єрiв для потреб будiвель-
I. М. Горбань 33
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
ної галузi. Через свої розмiри вони можуть iстотно
впливати як на характеристики рiчкової течiї, так
i на морфологiю русла, що може призвести до по-
шкоджень та руйнувань рiчкових iнженерних спо-
руд.
a
b
c
Рис. 11. Динамiка розмиву заглиблень рiзних
розмiрiв при w = 5 м, q = 10м2/c:
a – глибина = 1 м, дiаметр = 200 м;
b – глибина = 1 м, дiаметр = 400 м;
c – глибина = 2 м, дiаметр = 200 м
Тут аналiзується вплив на динамiку заглиблень
їхнiх розмiрiв. Розглядаються заглиблення пара-
болiчної форми глибиною 1 i 2 м та дiаметром 200
i 400 м. Довжина розрахункового каналу дорiв-
нює 2000 м, рiвень i витрата води в каналi скла-
дають w = 5 м та q = 10 м2/c вiдповiдно. Об’єм
стоку наносiв оцiнюється за степеневим законом
(7) при Ag = 0.005 с2/м, pm = 0.4. Параметри
розрахункової схеми: ∆x = 2.5 м, CFL = 0.5.
Результати розрахункiв (рис. 11) показують, що
якiсно динамiка розмиву пiщаного кар’єру не за-
лежить вiд його розмiрiв. Вiдносно транспорту на-
носiв у ньому є двi рiзнi областi: зона вiдкладен-
ня наносiв (верховий укiс) i зона розмиву (низовий
укiс). Пiдходячи до заглиблення, доннi наноси, що
транспортуються потоком, скочуються по верхо-
вому укосу. Частина з них вiдкладається на бiчнiй
стiнцi, через що її нахил зменшується. Пiсля того,
як ця стiнка набуває вертикальної форми, наноси
опускаються на дно, швидко засипаючи кар’єр. В
той самий час, удар потоку води об низовий укiс
викликає його розмив, внаслiдок чого збiльшує-
ться дiаметр кар’єру.
Крiм того, що в процесi розмиву змiнюється по-
перечний профiль заглиблення, воно ще й рухає-
ться вниз за течiєю. Порiвняльний аналiз резуль-
татiв, представлених на рис. 11, свiдчить про те,
що швидкiсть перемiщення заглиблень бiльше за-
лежить вiд їхньої глибини, нiж вiд поперечного
розмiру. Глибокi ями (рис. 11,c) рухаються наба-
гато повiльнiше, нiж мiлкi (рис. 11,a,b).
Рис. 12. Вiдносна помилка морфологiчних
розрахункiв при обтiканнi заглиблень рiзних розмiрiв:
1 – глибина = 1 м, дiаметр = 200 м; 2 – глибина = 1 м,
дiаметр = 400 м; 3 – глибина = 2 м, дiаметр = 200 м
Для оцiнки точностi проведених розрахункiв на
рис. 12 аналiзується вiдносна помилка (50), за
якою оцiнюються втрати твердої речовини в про-
цесi моделювання. Тут можна бачити, що у всiх
розглянутих випадках її порядок не перевищує
34 I. М. Горбань
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
10−3. Найбiльша помилка має мiсце при моделю-
ваннi розмиву глибоких ям, що пов’язане з вели-
кими градiєнтами донної поверхнi. Збiльшення то-
чностi розрахункiв може бути досягнуте зменше-
нням кроку просторової дискретизацiї.
Представленi результати свiдчать про високу
точнiсть та стiйкiсть виконаних гiдроморфоди-
намiчних розрахункiв, що пiдтверджує здатнiсть
розробленого чисельного методу моделювати ди-
намiку рiчкових течiй i транспорт наносiв у них на
великих часових перiодах. Очевидно, що успiх мо-
делювання реальних морфологiчних процесiв iсто-
тно залежить вiд гiпотез, покладених в основу мо-
делi, зокрема, про об’єм стоку наносiв та тертя
донної поверхнi. Степенева залежнiсть стоку на-
носiв вiд швидкостi потоку води, яка використо-
вувалася у вище наведених розрахунках, описує
процес лише наближено, i застосовується, голов-
ним чином, для тестування чисельних алгоритмiв
та якiсних оцiнок. Для одержання правдивих кiль-
кiсних характеристик морфологiчних процесiв не-
обхiдне залучення бiльш практичних оцiнок пото-
ку наносiв i гiдравлiчних параметрiв.
Рис. 13. Динамiка розмиву заглиблення глибиною 1 м
i дiаметром 200 м при w = 5 м, q = 5 м2/с
На рис. 13 наведенi результати моделювання
розмиву пiщаного кар’єру, одержанi з урахуван-
ням тертя на доннiй поверхнi, яке оцiнювалося за
законом Маннiнга (5). Значення коефiцiєнту шорс-
ткостi n бралося з дiапазону 0.0125÷ 0.025 с/м1/3,
який вiдповiдає даним натурних експериментiв на
руслових дiлянках рiчки Днiпро [32]. Потужнiсть
стоку наносiв оцiнювалася за формулою Мейера-
Петера i Мюллера (13), що добре описує об’єм сто-
ку наносiв для рiвнинних рiчок. Нахил дна в роз-
рахунках покладався 2◦. Дiаметр частинок пiску
ds = 0.4 мм, вiдношення густини седименту до гу-
стини води ρs/ρ = 2.65, що властиве для квар-
цу. Параметри гiдродинамiчної течiї: w = 5 м,
q = 5 м2/с. Отриманi результати дозволяють оцi-
нити характеристики "дрейфу" донного заглибле-
ння заданих розмiрiв (глибина = 1 м, дiаметр =
200 м). З рис. 13 випливає, що за добу таке за-
глиблення зноситься течiєю, середня швидкiсть
якої дорiвнює 1 м/c, приблизно на 20 м. Розви-
ненi моделi i виконанi розрахунки дають можли-
вiсть спрогнозувати розвиток гiдроморфодинамi-
чної ситуацiї в рiчкових басейнах за наявностi там
пiщаних кар’єрiв великих просторових масштабiв.
ВИСНОВКИ
1. Розвинений чисельний алгоритм для моде-
лювання гiдроморфодинамiчних процесiв у
вiдкритих водоймах. Вiн поєднує централь-
ну схему другого порядку Курганова-Ноелля-
Петрової для iнтегрування гiдродинамiчних
рiвнянь зi схемою "проти течiї" п’ятого по-
рядку, за якою розраховується еволюцiя дон-
ної поверхнi. Остання ґрунтується на зваже-
нiй iстотно неосцiлюючiй апроксимацiї потоку
наносiв.
2. Розглянутi одновимiрнi реалiзацiї цього ал-
горитму для незв’язаної i частково зв’язаної
систем гiдроморфодинамiчних рiвнянь. Ви-
конанi тестовi розрахунки еволюцiї пiщаного
пагорбу, розташованого на днi прямокутного
каналу. Порiвняння одержаних результатiв з
вiдомими аналiтичними розв’язками показу-
ють, що обидвi моделi забезпечують прийня-
тну точнiсть, але частково зв’язане моделюва-
ння краще пiдходить для довгострокових роз-
рахункiв, оскiльки майже повнiстю виключає
хибнi коливання шуканих розв’язкiв.
3. Проведенi довгостроковi розрахунки, в яких
отриманi часовi i просторовi характеристики
еволюцiї глибоких ям (кар’єрiв) та пiщаних
пагорбiв на рiчковому днi. Показано, що в
процесi розмиву не лише змiнюється попере-
чний профiль нерiвностi, а й вiдбувається її
перенос за течiєю.
4. Одержанi данi про швидкiсть переносу пiща-
них кар’єрiв у реальних умовах рiчкової течiї.
Вони дозволяють спрогнозувати розвиток гi-
дроморфодинамiчної ситуацiї в рiчкових ба-
сейнах зi складним дном, а також оцiнити
вплив крупномасштабних донних утворень на
еволюцiю течiй, розмивнi процеси поблизу бе-
регiв, стiйкiсть гiдротехнiчних споруд, тощо.
1. Сухомел Г. И. Исследование гидравлики
открытых русел и сооружений.– Киев: Нау-
кова думка, 1965.– 112 с.
I. М. Горбань 35
ISSN 1561 -9087 Прикладна гiдромеханiка. 2015. Том 17, N 1. С. 21 – 36
2. Гришанин К.В. Динамика русловых потоков.– Ле-
нинград: Госметеоиздат, 1979.– 276 с.
3. Wu W. Computational river dynamics.– London:
Taylor & Francis, 2008.– 488 p.
4. Hudson J., Sweby P.K. Formulations for numerically
approximating hyperbolic systems governing sedi-
ment transport // J. Scient. Comput.– 19.– 2003.–
P. 225–252.
5. Cunge J.A., Holly F.M., Verway A. Practical aspects
of computational river hydraulics.– London: Pitman,
1980.– 300 с.
6. Castro-Diaz M.J., Fernandez-Nieto E.D., Ferreiro
A.M. Sediment transport models in shallow water
equations and numerical approach by high order fini-
te volume methods // J. Comput. and Fluids.– 37.–
2008.– P. 299–316.
7. Goutiere L., Soares-Frazao S., Savary C., Laraichi T.,
Zech Y. One-dimensional model for transient flows
involving bed-load sediment transport and changes
in flow regimes // J. Hydr. Eng.– 134(6).– 2008.–
P. 726–735.
8. Kubatko E., Westerink J., Dawson C. An
unstructured grid morphodynamic model with
a discontinuous galerkin method for bed evolution //
J. Ocean Model.– 15.– 2006.– P. 71–89.
9. Bouchut F., Morales de Luna T. An entropy satisfyi-
ng scheme for two-layer shallow water equations with
uncoupled treatment // J. Math. Model. and Num.
Anal.– 42.– 2008.– P. 683–698.
10. Kurganov A., Noelle S., Petrova S. Semidiscrete
central-upwind schemes for hyperbolic conservation
laws and Hamilton-Jacobi equations // SIAM J. Sci.
Comput.– 23(3).– 2001.– P. 707–740.
11. Kurganov A., Lin C.-T. On the reduction of
numerical dissipation in central-upwind schemes //
Commun. Comput. Phys.– 2(1).– 2007.– P. 141–163.
12. Long W., Kirby J.T., Shao Z. A numerical scheme
for morphological bed level calculations // J. Coast.
Eng.– 55.– 2008.– P. 167–180.
13. de Saint-Venant B. Theorie du movement non-
permanent des eaux, avec application aux crues des
riviere at al’introduction des marees dans leur lit //
C. R. Acad. Sci., Paris.– 73.– 1871.– P. 147–154.
14. Einstein H.A. Formulas for the transportation of bed
load // J.Trans. ASCE.– 107.– 1942.– P. 561–573.
15. Grass A.J. Sediment transport by waves and
currents.– SERC London Center Marine Technology:
Report FL29, 1981.– 125 p.
16. Meyer-Peter E., Muller R. Formulas for bed-load
transport // Report on 2nd meeting on internati-
onal association on hydraulic structures research.–
Stockholm.– 1948.– P. 39–64.
17. Van Rijn Sediment transport (I): bed-load
transport // J. Hydroulic Div. ASCE.– 110.–
1984.– P. 1431–1456.
18. Nielsen P. Coastal bottom boundary layers and sedi-
ment transport.– London: World Scientific Pub Co
Ink, 1992.– 250 p.
19. Kurganov A., Tadmor E. New high-resolution
central schemes for nonlinear conservation laws and
convection-diffusion equations // J. Comput. Phys.–
160.– 2000.– P. 241–282.
20. Годунов С.К. Разностный метод для числен-
ного решения уравнений гидродинамики с ра-
зрывами // Математический сборник.– М: Изд-во
МГУ.– 1959.– С. 271–300.
21. Горбань В.О., Горбань I.М. Чисельнi моделi не-
стацiонарних руслових процесiв // Прикладна
гiдромеханiка.– 4.– 2013.– С. 55–75.
22. LeVeque R.J. Finite volume methods for hyperbolic
problems.– Cambridge: Cambridge University Press,
2002.– 130 p.
23. Lax P.D. Weak solutions of nonlinear hyperbolic
equations and their numerical computation //
Commun. Pur. Appl. Math.– 7.– 1954.– P. 159–193.
24. Sweby P. K. High resolution schemes using flux li-
miters for hyperbolic conservation laws // SIAM J.
Numer. Anal.– 21.– 1984.– P. 995–1011.
25. Louaked M., Hanich L. TVD scheme for the shallow
water equations // J. Hydraul. Res.– 27.– 1989.–
P. 321–332.
26. Gottlieb S., Shu C.W., Tadmor E. Strong stability-
preserving high order time discretization methods //
SIAM Rev.– 43.– 2001.– P. 89–112..
27. Kurganov A., Levy D. Central-upwind schemes for
the Saint-Venant system // ESAIM: M2AN.– 36.–
2002.– P. 397–425.
28. Johnson H.K., Zyserman J.A. Controlling spatial
oscillations in bed level update schemes // J. Coast.
Eng.– 46.– 2002.– P. 109–126.
29. Jensen J.H., Madsen E.Ø., Fredsøe J. Oblique flow
over dredged channels:II. Sediment transport and
morphology // J. Hydraul. Eng.– 125 (11).– 1999.–
P. 1190–1198.
30. Liu X.-D., Osher S., Chan T. Weighted essentially
non-oscillatory schemes // J. Comput. Phys.– 115.–
1994.– P. 200–211.
31. Jiang G.-S., Shu C.-W. Efficient implementation of
weighted ENO schemes // J. Comput. Phys.– 126.–
1996.– P. 202–214.
32. Исследование изменения режима р. Днепр и пер-
спективных методов добычи и транспортировки
песка в Каневском водохранилище.– Киев: Отчет
ИГМ АН УССР, 1988.– 250 с.
36 I. М. Горбань
|