Моделювання динаміки піщаних кар'єрів у руслових потоках

Розвинено двовимірний чисельний метод для розрахунку спільної еволюції водної течії та донної поверхні у відкритих водоймах, який враховує взаємний вплив гідродинамічного й морфологічного процесів. Він поєднує центральну схему "проти течії" Курганова-Ноелля-Петрової для інтегрування рівнян...

Повний опис

Збережено в:

Бібліографічні деталі
Дата:	2018
Автор:	Горбань, І.М.
Формат:	Стаття
Мова:	Ukrainian
Опубліковано:	Інститут гідромеханіки НАН України 2018
Назва видання:	Гідродинаміка і акустика
Онлайн доступ:	http://dspace.nbuv.gov.ua/handle/123456789/174286
Теги:	Додати тег Немає тегів, Будьте першим, хто поставить тег для цього запису!
Назва журналу:	Digital Library of Periodicals of National Academy of Sciences of Ukraine
Цитувати:	Моделювання динаміки піщаних кар'єрів у руслових потоках / І.М. Горбань // Гідродинаміка і акустика. — 2018. — Т. 1, № 2. — С. 132-159. — Бібліогр.: 24 назв. — укр.

Репозитарії

Digital Library of Periodicals of National Academy of Sciences of Ukraine

id	irk-123456789-174286
record_format	dspace
spelling	irk-123456789-1742862021-01-12T01:26:35Z Моделювання динаміки піщаних кар'єрів у руслових потоках Горбань, І.М. Розвинено двовимірний чисельний метод для розрахунку спільної еволюції водної течії та донної поверхні у відкритих водоймах, який враховує взаємний вплив гідродинамічного й морфологічного процесів. Він поєднує центральну схему "проти течії" Курганова-Ноелля-Петрової для інтегрування рівнянь мілкої води та зважену істотно неосцилюючу схему (WENO) п'ятого порядку, за якою розраховується еволюція донної поверхні. Потужність придонного потоку наносів оцінюється за степеневою формулою Грасса. Моделювання тестової задачі про розмив течією води конічного піщаного пагорба виявило довготермінову стійкість отриманого розв'язку, а також його відповідність даним інших дослідників. Развит численный метод для расчета совместной эволюции водного течения и донной поверхности в открытых водоемах, учитывающий взаимное влияние гидродинамического и морфологического процессов. Он сочетает центральную схему "против течения" Курганова-Ноелля-Петровой для интегрирования уравнений мелкой воды и взвешенную существенно неосциллирующую (WENO) схему пятого порядка, по которой рассчитывается эволюция донной поверхности. Мощность придонного потока наносов оценивается по степенной формуле Грасса. Моделирование тестовой задачи о размыве потоком воды конического песчаного холма выявило долговременную устойчивость полученного решения, а также его соответствие данным других авторов. A two-dimensional numerical method for simulating of coupled evolution of the water flow and the bed in open reservoirs is developed with the allowance for mutual influence of the morphological processes. It combines the Kurganov-Noelle-Petrova central upwind scheme for integration of shallow water equations and the fifth-order Weighted Essentially Non-Oscillatory (WENO) Scheme to describe the evolution of bed surface. The bedload sediment transport rate is estimated by power-law Grass formula. Simulation of the test problem on erosion of a conical sand dune in water flow revealed the longterm stability of the obtained solution, as well as its relevance to the data presented by other researchers. 2018 Article Моделювання динаміки піщаних кар'єрів у руслових потоках / І.М. Горбань // Гідродинаміка і акустика. — 2018. — Т. 1, № 2. — С. 132-159. — Бібліогр.: 24 назв. — укр. 2616-6135 DOI: doi.org/10.15407/jha2018.02.132 http://dspace.nbuv.gov.ua/handle/123456789/174286 532.517:627.157 uk Гідродинаміка і акустика Інститут гідромеханіки НАН України
institution	Digital Library of Periodicals of National Academy of Sciences of Ukraine
collection	DSpace DC
language	Ukrainian
description	Розвинено двовимірний чисельний метод для розрахунку спільної еволюції водної течії та донної поверхні у відкритих водоймах, який враховує взаємний вплив гідродинамічного й морфологічного процесів. Він поєднує центральну схему "проти течії" Курганова-Ноелля-Петрової для інтегрування рівнянь мілкої води та зважену істотно неосцилюючу схему (WENO) п'ятого порядку, за якою розраховується еволюція донної поверхні. Потужність придонного потоку наносів оцінюється за степеневою формулою Грасса. Моделювання тестової задачі про розмив течією води конічного піщаного пагорба виявило довготермінову стійкість отриманого розв'язку, а також його відповідність даним інших дослідників.
format	Article
author	Горбань, І.М.
spellingShingle	Горбань, І.М. Моделювання динаміки піщаних кар'єрів у руслових потоках Гідродинаміка і акустика
author_facet	Горбань, І.М.
author_sort	Горбань, І.М.
title	Моделювання динаміки піщаних кар'єрів у руслових потоках
title_short	Моделювання динаміки піщаних кар'єрів у руслових потоках
title_full	Моделювання динаміки піщаних кар'єрів у руслових потоках
title_fullStr	Моделювання динаміки піщаних кар'єрів у руслових потоках
title_full_unstemmed	Моделювання динаміки піщаних кар'єрів у руслових потоках
title_sort	моделювання динаміки піщаних кар'єрів у руслових потоках
publisher	Інститут гідромеханіки НАН України
publishDate	2018
url	http://dspace.nbuv.gov.ua/handle/123456789/174286
citation_txt	Моделювання динаміки піщаних кар'єрів у руслових потоках / І.М. Горбань // Гідродинаміка і акустика. — 2018. — Т. 1, № 2. — С. 132-159. — Бібліогр.: 24 назв. — укр.
series	Гідродинаміка і акустика
work_keys_str_mv	AT gorbanʹím modelûvannâdinamíkipíŝanihkarêrívuruslovihpotokah
first_indexed	2025-07-15T11:13:30Z
last_indexed	2025-07-15T11:13:30Z
_version_	1837711240841920512
fulltext	ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. УДК 532.517:627.157 МОДЕЛЮВАННЯ ДИНАМIКИ ПIЩАНИХ КАР’ЄРIВ У РУСЛОВИХ ПОТОКАХ I. М. Горбань† Iнститут гiдромеханiки НАН України вул. Желябова, 8/4, 03057, Київ, Україна †E-mail: ivgorban@gmail.com Получено 09.11.2017 Розвинено двовимiрний чисельний метод для розрахунку спiльної еволюцiї водної течiї та донної поверхнi у вiдкритих водоймах, який враховує взаємний вплив гiд- родинамiчного й морфологiчного процесiв. Вiн поєднує центральну схему «проти течiї» Курганова—Ноелля—Петрової для iнтегрування рiвнянь мiлкої води та зва- жену iстотно неосцилюючу схему (WENO) п’ятого порядку, за якою розраховуєть- ся еволюц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д 20∘ до 30∘. Отриман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щаний кар’єр, кут розмиву 132 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. 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чних споруд. Морфолог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ї [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лкої води для опису течiї та рiвняння збереження маси наносiв, за яким розраховується еволюцiя дна. Для її замикання використовуються емпiричнi залежностi, якi пов’язують потужнiсть потоку наносiв з характеристиками водної течiї. Докладний огляд цих законiв представлено в роботах [22,33]. При моделюваннi значних деформацiй донної поверхнi, якi цiкавитимуть нас у цьо- му дослiдженнi, необхiдно враховувати взаємний вплив водної течiї й донних процесiв, а це потребує застосування зв’язаного пiдходу, коли гiдродинамiчнi й морфологiчнi рiв- няння розв’язуються одночасно. В роботi [44] показано, що при зв’язаному 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д розвивається до- сить активно [22,55], одержанi розв’язки стосуються одновимiрних моделей зi степеневою залежнiстю потоку наносiв вiд швидкостi водної течiї. У цьому дослiдженнi застосовується стратегiя розщеплення гiдроморфодинамiчної системи за часом, коли рiвняння мiлкої води розв’язуються при фiксованiй формi дна, пiсля чого вона оновлюється на тому ж часовому кроцi. Цей алгоритм простiший, нiж 133 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. повн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внянь [44, 66, 77]. Оскiльки результати натурних дослiджень пiщаних кар’єрiв вказують на їхню складну i довготривалу еволюцiйну iсторiю [88], то для моделювання таких утворень слiд засто- совувати чисельнi пiдходи, якi здатнi забезпечити стiйкiсть i точнiсть морфологiчних розрахункiв протягом довгого часу. В роботi [77] зроблено пор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вняння — WENO-approximation. WENO-схеми вперше були представленi в роботi [99] i з того часу активно залучаються до чисельного iнтегрування диференцiйних рiвнянь у частинних похiдних, зокрема, рiвнянь гiперболiчного типу. Їхнi переваги полягають у високому порядку точностi в областях монотонностi шуканих величин i стiйкiй неос- цилюючiй поведiнцi розв’язкiв поблизу розривiв. У роботi [33] WENO-схему було 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 застосування WENO-апроксимацiї. Як i в роботi [33], для просторової дис- кретизацiї рiвнянь мiлкої води тут використано прямокутну сiтку й центральну схему з числовим потоком у формi Курганова—Ноелля—Петрової [1010] та антидифузiйним чле- ном у формi Кургановi—Лiня [1111]. Еволюцiя гiдродинамiчного процесу в часовi розра- ховувалась методом Рунге—Кутта третього порядку, який належить до SSP-схем, що забезпечують стiйкий розв’язок [1212]. Морфологiчне рiвняння iнтегрувалося за часом методом Ейлера, що значно спрощує застосування моделi до розв’язання практичних проблем. Розвинений алгоритм було протестовано на вiдомiй модельнiй задачi про розмив по- током цилiндричного пiщаного пагорба, розташованого на днi прямокутного каналу [1313]. При цьому виявлено стiйкiсть отримуваних розв’язкiв на довгому часовому промiжку, а також їхню вiдповiднiсть даним iнших авторiв [1313, 1414]. Це свiдчить про придатнiсть запропонованої чисельної схеми для розв’язання складних фiзично обґрунтованих мор- фодинамiчних проблем. 134 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. Проведено розрахунки довгострокових часових i просторових характеристик ево- люцiї пiщаних кар’єрiв на рiчковому днi. Результати моделювання продемонстрували утворення складного рельєфу на доннiй поверхнi та його пересування вниз за течiєю. Так, цилiндричне заглиблення, яке на початку мало локалiзований характер, у процесi еволюцiї розпадалось на систему ям i невеликих насипiв, якi iстотно поширювались уздовж поперечного перерiзу каналу. Отриманi величини кута розмиву, який характе- ризує поширення заглиблення на днi, знаходиться в межах вiд 20∘ до 30∘ 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чок тощо. 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лкої води набуде такого вигляду [11]: ℎ𝑡 + (𝑢ℎ)𝑥 + (𝑣ℎ)𝑦 = 0, (𝑢ℎ)𝑡 + [︂ ℎ𝑢2 + 1 2 𝑔ℎ2 ]︂ 𝑥 + (ℎ𝑢𝑣)𝑦 = −𝑔ℎ(𝑆0𝑥 − 𝑆𝑓𝑥), (𝑣ℎ)𝑡 + (ℎ𝑢𝑣)𝑥 + [︂ ℎ𝑣2 + 1 2 𝑔ℎ2 ]︂ 𝑦 = −𝑔ℎ(𝑆0𝑦 − 𝑆𝑓𝑦), (1) де 𝑡 — час; 𝑥 i 𝑦 — поздовжня й поперечна координати вiдповiдно; 𝑔 — прискорення вiльного падiння; ℎ — глибина води; 𝑢i 𝑣 — компоненти швидкостi потоку в 𝑥- i 𝑦- напрямках вiдповiдно; 𝑆0 = (𝑆0𝑥, 𝑆0𝑦) — геометричний нахил донної поверхнi; 𝑆𝑓 = (𝑆𝑓𝑥, 𝑆𝑓𝑦) — нахил тертя. Рiвняння для визначення форми донної поверхнi випливає з закону збереження маси твердої речовини в рiчковому потоцi. Для придонного потоку наносiв воно вiдоме як рiвняння Екснера [22]: (1 − 𝑝𝑚) 𝜕𝑧𝑏 𝜕𝑡 + 𝜕𝑞1 𝜕𝑥 + 𝜕𝑞2 𝜕𝑦 = 0, (2) де 𝑝𝑚 — пористiсть донного ґрунту; 𝑧𝑏 = 𝑧𝑏(𝑥, 𝑦, 𝑡) — функцiя дна; 𝑞1, 𝑞2 — потужностi потоку наносiв в 𝑥- i 𝑦-напрямках вiдповiдно. Зазначимо, що рiвняння (2)(2) записане в припущеннi однорiдностi дна, коли пористiсть поверхнi й розмiри седиментних частинок скрiзь однаковi. 135 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. а б Рис. 1. Схема руслового потоку з розмивним дном: а — в площинi 𝑂𝑥𝑧, б — в площинi 𝑂𝑥𝑦 Параметри системи рiвнянь (1)(1), (2)(2) показанi на Рис. 1Рис. 1, який схематично iлюструє руслову течiю в площинах 𝑂𝑥𝑧 i 𝑂𝑥𝑦. Невiдомими тут є глибина води ℎ, витрати води 𝑄𝑥 = 𝑢ℎ i 𝑄𝑦 = 𝑣ℎ, а також функцiя дна 𝑧𝑏. Геометричний нахил поверхнi становить собою градiєнт функцiї 𝑧𝑏, так що 𝑆0𝑥 = −𝜕𝑧𝑏 𝜕𝑥 , 𝑆0𝑦 = −𝜕𝑧𝑏 𝜕𝑦 . (3) Нахил тертя розраховується за емпiричним законом Маннiнга 𝑆𝑓𝑥 = 𝑛2𝑢 √ 𝑢2 + 𝑣2 ℎ4/3 , 𝑆𝑓𝑦 = 𝑛2𝑣 √ 𝑢2 + 𝑣2 ℎ4/3 , (4) де 𝑛 — коефiцiєнт шорсткостi донної поверхнi, значення якого для конкретних водойм визначається в натурних експериментах. Величини 𝑞1, 𝑞2 залежать вiд типу твердої речовини й характеристик гiдродинамiч- ної течiї. Вiдомi оцiнки для їх визначення сформульовано здебiльшого для гранульова- них частинок, i вони мають вигляд ймовiрнiсних або детермiнiстичних законiв [22, 33]. У цьому дослiдженнi обмежимося застосуванням найпростiшої степеневої формули Грас- са, в який потужнiсть потоку наносiв виражається безпосередньо через швидкiсть вод- ної течiї: 𝑞1 = 𝐴𝑔𝑢(𝑢2 + 𝑣2)(𝑚−1)/2, 𝑞2 = 𝐴𝑔𝑣(𝑢2 + 𝑣2)(𝑚−1)/2. (5) Тут 𝐴𝑔 i 1 ≤ 𝑚 ≤ 4 — емпiричнi константи. Вибiр їхнiх величин залежить вiд розмiру частинок, що переносяться, в’язкостi рiдини та iнших фiзичних факторiв. Значення 𝐴𝑔 ≪ 1 вiдповiдають слабкiй взаємодiї мiж течiєю води й транспортом наносiв. Якщо ж 𝐴𝑔 ∼ 1, то має мiсце сильна взаємодiя. 3. ЧИСЕЛЬНА СХЕМА Для iнтегрування системи рiвнянь (1)(1), (2)(2) застосуємо стратегiю розщеплення гiд- роморфодинамiчного процесу в рамках часового промiжку ∆𝑡, коли рiвняння мiлкої води розв’язуються при фiксований формi дна, а розрахованi значення гiдродинамiч- них змiнних використовуються у морфологiчному рiвняннi, щоб оновити топографiю 136 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. донної поверхнi. Протягом останнього часу цей пiдхiд активно розвивається [66], оскiль- ки вiн значно простiший, нiж повнiстю зв’язаний i, водночас, дозволяє врахувати вза- ємний вплив гiдродинамiчного й морфологiчного процесiв навiть при швидких змiнах параметрiв. Розроблений чисельний алгоритм поєднує центральну схему другого порядку Кур- ганова—Ноелля—Петрової (KNP-scheme) [1010,1111], використану для iнтегрування рiвнянь мiлкої води, зi схемою «проти течiї» п’ятого порядку, за якою розраховується еволюцiя донної поверхнi. Остання ґрунтується на зваженiй iстотно неосцилюючiй апроксима- цiї (WENO-approximation) потоку наносiв, запропонованiй в роботi [77]. Одновимiрна реалiзацiя описаного алгоритму представлена в роботi [33], де на багатьох тестових при- кладах показано, що вiн повнiстю виключає нефiзичнi осциляцiї шуканих розв’язкiв, завдяки чому дозволяє успiшно моделювати довготривалi морфологiчнi процеси. Адап- тацiя розвиненого чисельного методу до двовимiрного випадку значно розширює межi його застосування, зокрема, дозволяє моделювати еволюцiю складних форм донної по- верхнi. 3.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й, а й зберiгати їх протягом довгого часу. Од- ним iз найпростiших i найефективнiших пiдходiв, якi вiдповiдають згаданим вимогам, є центральна схема «проти течiї», обґрунтована в роботах [1010,1111,1515]. Вона належить до проекцiйно-еволюцiйних методiв типу Годунова, якi мають високу роздiльну здатнiсть i дозволяють знаходити розривнi розв’язки. Докладне описання напiвдискретної центральної схеми «проти течiї» представлене в публiкацiях [1010, 1616], що дозволяє обмежитись перелiком лише загальних її особливо- стей. Так, для просторової дискретизацiї рiвнянь мiлкої води застосовуються прямо- кутнi скiнченнi об’єми. Виходячи з цього, замiсть локальних фiзичних величин роз- глядаються їхнi значення, усередненi по елементарному об’єму. Реконструкцiя змiнних на межах комiрки виконується за допомогою кусково-лiнiйної апроксимацiї. Для ро- зрахунку похiдних використовується minmod-лiмiтер, що дозволяє усунути осциляцiї шуканих величин в областях з великими градiєнтами [1717]. Розриви змiнних, що вно- сяться в схему в процесi реконструкцiї, оцiнюються на основi iнформацiї про локальнi хвильовi швидкостi — власнi числа якобiана дослiджуваної системи. Конвективнi члени у лiвих частинах двох останнiх рiвнянь системи (1)(1) апроксиму- ються дискретними операторами (чисельними потоками), якi розраховуються на основi реконструйованих змiнних на границях елементiв сiтки. У цьму дослiдженнi викори- стовується чисельний потiк у формi Курганова—Ноелля—Петрової [1010], побудований таким чином, що iнтегрування за часом виконується окремо по областi розриву i по тiй частинi об’єму, де змiннi неперервнi. Ширина розриву обчислюється за допомогою одностороннiх мiнiмальної та максимальної хвильових швидкостей, завдяки чому за- 137 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. пропоновану схему, хоч вона й залишається центральною за природою, слiд вiднести до схем проти течiї з низькою чисельною дифузiєю. Для подальшого зменшення си- стемної дифузiї в роботi [1111] запропоновано ввести антидифузiйний член, за допомогою якого забезпечується бiльш якiсне виконання проекцiйної процедури, в ходi якої проiн- тегрованi за часом реконструйованi змiннi перераховуються на початкову сiтку, тобто усереднюються по комiрцi. Зазначимо, що використана просторова апроксимацiя має другий порядок точностi. Одержанi шляхом просторової дискретизацiї звичайнi диференцiйнi рiвняння iнте- груються за часом за допомогою стiйкого солвера Рунге—Кутта третього порядоку точ- ностi (SSP-RK), вигляд якого можна знайти, наприклад, у роботi [1616]. Оскiльки цей алгоритм є явним за часом, його стiйкiсть повинна забезпечуватись виконанням умови Куранта—Фрiдрiхса—Левi: ∆𝑡 min{∆𝑥,∆𝑦} max{𝑎max, 𝑏max} ≤ CFL, (6) де 𝑎max i 𝑏max — максимальнi локальнi швидкостi поширення збурень у розрахунковiй областi в 𝑥- та 𝑦-напрямках вiдповiдно; CFL — число Куранта—Фрiдрiхса—Левi, яке є мiрою поширення збурень за часовий промiжок ∆𝑡. У роботi [1010] показано, що стiйкiсть алгоритму забезпечується при CFL ≤ 0.5. Застосування KNP-схеми до неоднорiдних законiв збереження, до яких належать рiвняння мiлкої води, представлене в [1818]. Консервативнiсть, яка потребує збережен- ня стацiонарних розв’язкiв з порядком точностi не нижче другого, забезпечується тут застосуванням кусково-лiнiйної апроксимацiї дна i спецiальних квадратур для дискре- тизацiї правої частини рiвнянь (1)(1). Збереження позитивної глибини (positive preserving property) у всiй областi досягається завдяки переходу в рiвняннях вiд змiнної ℎ до змiн- ної 𝑤, яка описує форму вiльної поверхнi, а також за допомогою спецiальної процедури корекцiї розв’язкiв. Виконанi в статтi [1616] докладн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 дiлянок незмоченної донної поверхнi. 3.2. Дискретизацiя морфологiчного рiвняння Для кращого розумiння природи еволюцiйного рiвняння (2)(2) запишемо його у виглядi 𝜕𝑧𝑏 𝜕𝑡 + 𝜉 (︂ 𝜕𝑞1 𝜕𝑥 + 𝜕𝑞2 𝜕𝑦 )︂ = 0, (7) 138 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. де 𝜉 = 1/(1 − 𝑝𝑚). Враховуючи, що поле течiї, а значить, i потужностi потоку наносiв залежать вiд форми донної поверхнi, рiвняння (7)(7) може розглядатися як нелiнiйний закон збережен- ня вiдносно змiнної 𝑧𝑏. Тому, для його розв’язання повиннi застосовуватись методи, якi можуть знаходити розривнi розв’язки (shock capturing schemes). Вони ґрунтуються на хвильовiй апроксимацiї еволюцiйного рiвняння, яка випливає з припущення про те, що потужностi придонного потоку наносiв 𝑞1, 𝑞2 залежать лише вiд форми донної поверх- нi 𝑧𝑏. Його коректнiсть обумовлена тим, що просторова дискретизацiя морфологiчного рiвняння виконується в межах окремого кроку за часом, коли, за виключенням форми дна, всi фактори, якi можуть впливати на потiк наносiв — глибина i швидкiсть течiї, фiзичнi властивостi частинок, — фiксованi [1919]. Поклавши 𝜕𝑞1 𝜕𝑥 = 𝜕𝑞1 𝜕𝑧𝑏 𝜕𝑧𝑏 𝜕𝑥 , 𝜕𝑞2 𝜕𝑦 = 𝜕𝑞2 𝜕𝑧𝑏 𝜕𝑧𝑏 𝜕𝑦 , (8) отримуємо з (7)(7) рiвняння адвекцiї вiдносно змiнної 𝑧𝑏, згiдно з яким морфологiчна еволюцiя вiдбувається як нелiнiйне поширення деформацiй дна у напрямку переносу домiшок: 𝜕𝑧𝑏 𝜕𝑡 + 𝐶𝑥(𝑧𝑏) 𝜕𝑧𝑏 𝜕𝑥 + 𝐶𝑦(𝑧𝑏) 𝜕𝑧𝑏 𝜕𝑦 = 0. (9) Тут через 𝐶𝑥(𝑧𝑏) = 𝜉 𝜕𝑞1 𝜕𝑧𝑏 , 𝐶𝑦(𝑧𝑏) = 𝜉 𝜕𝑞2 𝜕𝑧𝑏 позначено фазовi швидкостi поширення збурень на доннiй поверхнi в поздовжньому й поперечному напрямках вiдповiдно. Методи дискретизацiї рiвняння адвекцiї за часом достатньо широко представленi в лiтературi — вiд найпростiшої схеми Ейлера до багатокрокових SSP-схем високих порядкiв, розвинутих на основi метода Рунге—Кутта [1212]. Головна складнiсть при побу- довi чисельних апроксимацiй для конвективно-дифузiйних проблем, так само, як i для гiперболiчних законiв збереження, пов’язана з просторовою дискретизацiєю, оскiльки при цьому мають враховуватись розриви шуканих величин. В роботi [77] показано, що дискретизацiя рiвняння (9)(9) класичними центральними алгоритмами типу Лакса—Вендорфа призводить до високочастотних осциляцiй уро- зв’язках навiть при застосуваннi регуляризацiйних процедур, таких як введення штуч- ної в’язкостi [1919] або фiльтрування [2020]. Схеми проти течiї, якi будуються за допомогою рiманiвських солверiв, виявляють бiльшу стiйкiсть, нiж центральнi, але мають низь- кий порядок точностi, тодi як акуратне моделювання морфологiчних процесiв потребує схем високого порядку. Сучаснi алгоритми просторової дискретизацiї диференцiйних операторiв гiперболiч- ного типу ґрунтуються на iстотно неосцилюючих схемах (ENO або WENO), якi забезпе- чують високiй порядок точностi в областях монотонностi шуканих величин i стiйкiсть розв’язкiв поблизу розривiв [99, 2121]. В роботi [77] зважену iстотно неосцилюючу схему було адаптовано до розв’язання одновимiрного морфологiчного рiвняння. У цьому до- слiдженнi її розвинуто для двовимiрного випадку, а також iнкорпоровано до загальної чисельної моделi, у рамках якої виконується зв’язане iнтегрування гiдродинамiчної й морфологiчної задач. 139 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. Суть цього методу полягає в тому, що кожна зi складових транспорту наносiв роз- кладається на два потоки — позитивний i негативний. Так, для величини 𝑞1 можемо записати 𝑞1 = 𝑞−1 + 𝑞+1 , (10) де 𝑞−1 i 𝑞+1 оцiнюються за допомогою одностороннiх фазових швидкостей: 𝐶− 𝑥 = min {𝐶𝑥, 0}, 𝐶+ 𝑥 = max {𝐶𝑥, 0}. (11) З (11)(11) випливає, що 𝐶𝑥 = 𝐶− 𝑥 + 𝐶+ 𝑥 . Тодi, для одностороннiх потокiв, що входять до (10)(10), маємо: 𝑞−1 = 1 𝜉 𝑧𝑏∫︁ 0 𝐶− 𝑥 (𝑧)𝑑𝑧, 𝑞+1 = 1 𝜉 𝑧𝑏∫︁ 0 . (12) Аналогiчно оцiнюється потужнiсть переносу домiшок у поперечному напрямку 𝑞2. Для дискретизацiї рiвняння (7)(7) за просто- Рис. 2. Елементарний об’єм розрахункової сiтки ром на розрахункову область накладається рiвномiрна сiтка, яка складається з елементар- них об’ємiв: Ω𝑖𝑗 = {𝜉, 𝜂 : \|𝜉 − 𝑥𝑖\| < ∆𝑥/2, \|𝜂 − 𝑦𝑖\| < ∆𝑦/2}, де 𝑖 = 1, 2, . . . , 𝑁𝑥; 𝑗 = 1, 2, . . . , 𝑁𝑦; 𝑁𝑥 i 𝑁𝑦 — кiлькiсть елементiв сiтки в поздовжньому i по- перечному напрямках вiдповiдно (Рис. 2Рис. 2). В чисельнiй схемi дно буде оновлюватися у вуз- лових точках (𝑥𝑖, 𝑦𝑗), а потоки твердої речови- ни 𝑞1 i 𝑞2 оцiнюватися на межах елементарних об’ємiв. При застосуваннi формули Ейлера для iн- тегрування процесу за часом чисельна схема для розрахунку форми донної поверхнi набуває такого вигляду: (𝑧𝑏) 𝑛+1 𝑖,𝑗 − (𝑧𝑏) 𝑛 𝑖,𝑗 ∆𝑡𝑚 = −𝜉 (︂ (𝑞1)𝑖+1/2,𝑗 − (𝑞1)𝑖−1/2,𝑗 ∆𝑥 + (𝑞2)𝑖,𝑗+1/2 − (𝑞2)𝑖,𝑗−1/2 ∆𝑦 )︂ . (13) Тут ∆𝑡𝑚 — крок дискретизацiї за часом морфологiчного процесу; величини (𝑞1)𝑖±1/2,𝑗 i (𝑞2)𝑖,𝑗±1/2 — чисельнi апроксимацiї складових транспорту наносiв 𝑞1, 𝑞2 на кiнцях еле- ментарного об’єму. Очевидно, що якiсть результатiв, обчислених за схемою (13)(13), залежить вiд того, наскiльки точно оцiнюються величини (𝑞1)𝑖±1/2,𝑗 i (𝑞2)𝑖,𝑗±1/2. В iстотно неосцилюючих алгоритмах вони розщеплюються на позитивний i негативний потоки, як показано на Рис. 2Рис. 2. У вiдповiдностi до формули (10)(10) маємо (𝑞1)𝑖±1/2,𝑗 = (𝑞−1 )𝑖±1/2,𝑗 + (𝑞+1 )𝑖±1/2,𝑗, (𝑞2)𝑖,𝑗±1/2 = (𝑞−2 )𝑖,𝑗±1/2 + (𝑞+2 )𝑖,𝑗±1/2. (14) 140 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. Одностороннi потоки в (14)(14) обчислюються, виходячи зi значень 𝑞1 i 𝑞2 в сусiднiх iз заданою межею вузлових точках сiтки. При цьому враховується знак фазової швидкостi на цiй межi, який визначає напрямок поширення збурень донної поверхнi. Докладний опис WENO-схеми для одновимiрного морфологiчного рiвняння можна знайти в роботах [33,77]. Її адаптацiя до двовимiрного випадку виконується автоматично, тому її деталi тут не наводяться. Зазначимо лише, що ця схема має п’ятий порядок точностi, а завдяки врахуванню одностороннiх хвильових швидкостей вона належить до схем «проти течiї», завдяки чому здатна забезпечити високу стiйкiсть розрахункiв. У статтi [77] показано, що застосування багатокрокових схем типу Рунге—Кутта до iнтегрування морфологiчного рiвняння за часом майже не впливає на результати, але водночас потребує перерахунку гiдродинамiчних параметрiв течiї на промiжних часо- вих кроках. Такий пiдхiд може виявитись занадто затратним при виконаннi практичних розрахункiв. Представлений тут алгоритм, який поєднує WENO-схему з методом Ейле- ра (WENO-ER), є оптимальним з точки зору досягнення балансу точностi результатiв i простоти реалiзацiї. Зазначимо також, оскiльки схема (13)(13) — явна за часом, вона потребує виконання стандартного критерiю стiйкостi Куранта—Фрiдрiхса—Левi: 𝑐∆𝑡𝑚 min{∆𝑥,∆𝑦} ≤ 1, 𝑐 = max(\|𝐶𝑥𝑖+1/2,𝑗\|, \|𝐶𝑦𝑖,𝑗+1/2\|). (15) Втiм, при зв’язаному iнтегруваннi гiдроморфодинамiчної системи (1)(1), (2)(2), для вибору кроку за часом необхiдно реалiстично оцiнити характернi часовi масштаби водного по- току й транспорту наносiв. Очевидно, що гiдродинамiчнi змiни розвиваються набагато швидше, нiж морфологiчнi, тому в нашiй моделi крок за часом обирається з критерiю (6)(6). При цьому нiяких умов на характер взаємодiї мiж водною течiєю i дном не накла- дається. Зазначимо також, що у всiх морфологiчних розрахунках перевiряється збереження об’єму твердої речовини. Загальний об’єм пiщаних наносiв в областi може бути вираху- ваний за формулою 𝑉 (𝑡) = +∞∫︁ −∞ +∞∫︁ −∞ (𝑧𝑏(𝑥, 𝑦, 𝑡) − 𝑧𝑏∞) 𝑑𝑥𝑑𝑦, (16) де 𝑧𝑏∞ — умовний рiвень незбуреного дна. Вiдносна похибка розрахункiв щодо вико- нання закону збереження речовини обчислюється як Error = 𝑉 (𝑡) − 𝑉 (0) 𝑉 (0) . (17) 4. ТЕСТОВI РОЗРАХУНКИ У цьому роздiлi на вiдомих тестових прикладах буде оцiнено здатнiсть представ- леної чисельної моделi коректно розраховувати гiдродинамiчнi й морфологiчнi змiни у вiдкритих водоймах. Спочатку розглянемо застосування кожної зi схем — KNP чи WENO-ER — до вiдповiдного процесу, а потiм проiлюструємо результати, отриманi при зв’язаному моделюваннi гiдроморфодинамiчної задачi. 141 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. 4.1. Малi збурення вiльної поверхнi над неплоским нерозмивним дном Цю тестову задачу запропоновано в роботi [2222] для апробацiї чисельних схем, за яки- ми iнтегруються рiвняння мiлкої води. У нiй оцiнюється здатнiсть алгоритму корект- но описувати хвилi, що генеруються нерiвним дном. Розрахункова область має форму прямокутного каналу розмiрами 2× 1 м з сильно вигнутою донною поверхнею (Рис. 3Рис. 3). Конфiгурацiя дна задається функцiєю 𝑧(𝑥, 𝑦) = 0.8𝑒[−5(𝑥−0.9)2−50(𝑦−0.5)2]. Рiвень води в каналi скрiзь дорiвнює 1 м, за винятком малого локалiзованого збу- Рис. 3. Конфiгурацiя донної поверхнi Рис. 4. Еволюцiя малої хвилi над вигнутим нерозмивним дном: а — 𝑡 = 0.24 с, б — 𝑡 = 0.36 с, в — 𝑡 = 0.48 с, г — 𝑡 = 0.60 с 142 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. рення висотою 1 см, фронт якого паралельний осi 𝑂𝑦: 𝑤(𝑥, 𝑦, 0) = {︃ 1.01 м, якщо 𝑥 ∈ [0.05, 0.15]; 1 м, якщо 𝑥 ∈ [0, 0.05) ∪ (0.15, 2]. Рiдина на початку процесу нерухома, тобто для складових швидкостi маємо нульовi початковi умови: 𝑢(𝑥, 𝑦, 0) = 0, 𝑣(𝑥, 𝑦, 0) = 0. На Рис. 4Рис. 4 представленi одержанi в розрахунках трансформацiї вiльної поверхнi про- тягом 0.6 с (∆𝑥 = ∆𝑦 = 0.01 м). З графiкiв видно , що внаслiдок початкового збурення поверхнi над пагорбом утворюється хвиля, яка рухається в додатному напрямку осi 𝑂𝑥. При її проходженнi над вершиною пагорба утворюється вирва, а сама вона роз- падається на систему вторинних хвиль, якi поступово затухають. Отриманi результати добре узгоджуються з даними iнших авторiв, зокрема, представленими в дисертацiйнiй роботi [2323], де для моделювання процесу застосовувалася схема «проти течiї». 4.2. Розмив пiщаного пагорба в стацiонарному потоцi Цей класичний приклад застосовується для апробацiї чисельних схем, якi моделюють транспорт наносiв в придоннiй областi [77]. Розрахуємо розмив симетричного пiщаного пагорба потоком води зi сталими параметрами. Початкова форма дна задається функ- цiєю 𝑧𝑏(𝑥, 0) = −ℎ0 + 2 exp [−𝛽(𝑥− 𝑥𝑐)] , де ℎ0 — глибина води над горизонтальним дном; 𝑥𝑐 — центр пагорба; 𝛽 — параметр, який визначає крутизну пагорба. Покладемо ℎ0 = 6 м, 𝑥𝑐 = 150 м, 𝛽 = 0.01. Довжина розрахункової областi дорiвнює 300 м. Рiвень та витрата води протягом всього розрахунку вважаються незмiнними: 𝑤(𝑥, 𝑡) = 0, 𝑞(𝑥, 𝑡) = 10 м2/с. Потужнiсть потоку зважених части- Рис. 5. Розмив пiщаного пагорба: 1 — 𝑡 = 0, 2 — 𝑡 = 600 c, 3 — 𝑡 = 2000 c, 4 — 𝑡 = 6000 c нок в розрахунках оцiнюється за степе- невою формулою (5)(5) при 𝐴𝑔 = 0.001, пористiсть ґрунту — 𝑝𝑚 = 0.4. В роз- глянутому випадку вибрано такi пара- метри розрахункової схеми: ∆𝑥 = 1 м, ∆𝑡 = 0.1 с. Результати розрахункiв та їхнє порiвняння з точним розв’язком [77] представленi на Рис. 5Рис. 5. З графiка видно, що розмив пiщаного пагорба супроводжується генерацiєю гострого фронту седиментних частинок, який рухається вниз за течiєю. Важливо, що чисельна схема описує змiни конфiгурацiї дна без хибних осциляцiй i не втрачає стiйкостi протягом достатньо довгого розрахункового перiоду. Отриманi результати добре узгоджуються з точним розв’язком. 143 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. 4.3. Розмив пiщаного пагорба в нестацiонарному потоцi Для верифiкацiї чисельної схеми, яка виконує зв’язане моделювання гiдродинамiчного й морфологiчного процесiв у двовимiрнiй постановцi, розглнемо тестову задача про розмив пiщаного пагорба [1313, 1414]. Тут моделюються як еволюцiя донної поверхнi, так i викликанi нею змiни гiдродинамiчних параметрiв течiї. Розрахункова область являє собою прямокутний канал з горизонтальними розмiрами 1500 × 1000 м, на днi якого розташований пiщаний пагорб у формi конуса висотою 1 м i дiаметром 200 м. Таким чином, форма дна задається функцiєю (всi розмiри в метрах): 𝑧𝑏(𝑥, 𝑦, 0) = ⎧⎨⎩ sin2 (︂ 𝜋(𝑥− 300) 200 )︂ sin2 (︂ 𝜋(𝑦 − 400) 200 )︂ при 𝑥 ∈ [300, 500], 𝑦 ∈ [400, 600]; 0 в iнших випадках. У початковий момент поверхня води незбурена, її рiвень становить 10 м, а витрата води — 10 м3/с. Виходячи з цього, початковi умови для гiдродинамiчних параметрiв мають вигляд ℎ(𝑥, 𝑦, 0) = 10 − 𝑧𝑏(𝑥, 𝑦, 0), 𝑢(𝑥, 𝑦, 0) = 10 ℎ(𝑥, 𝑦, 0) , 𝑣(𝑥, 𝑦, 0) = 0. На входi розрахункової областi задаються глибина i витрата води, на а виходi накла- дається умова вiдкритого потоку. На бiчних стiнках, паралельних основному потоку, накладається умова вiдбиття, згiдно з якою рiдина i наноси не перетiкають через них. Оскiльки розглядається докритична течiя (Fr = √ 𝑢2 + 𝑣2/(𝑔ℎ) ≈ 0.1), наведених умов достатньо для коректної постановки граничної задачi для системи рiвнянь (1)(1), (2)(2). До- кладний опис процедури iмплементацiї граничних умов в дискретну схему з кiнцевими об’ємами наведено в статтi [1616]. Розрахункова область дискретизується вiдрiзками ∆𝑥 = ∆𝑦 = 20 м, а крок за часом обчислюється з критерiю Куранта—Фрiдрiхса—Левi (6)(6) при 𝐶𝐹𝐿 = 0.5. Потужнiсть потоку седиментних частинок оцiнюється за степеневою формулою (5)(5) при 𝑚 = 3, а також 𝐴𝑔 = 0.001 або 1 (це що вiдповiдає слабкiй i сильнiй взаємодiї транспорту наносiв з водною течiєю). Пористiсть ґрунту становить 𝑝𝑚 = 0.4. Для уникнення рiзкого старту морфологiчного процесу на першому етапi моделюєть- ся лише гiдродинамiчна течiя до досягнення нею стацiонарного стану при заданiй формi дна. Останнiй визначається з умови \|𝑤𝑛+1 𝑖,𝑗 − 𝑤𝑛 𝑖,𝑗\| ≤ 𝜀 ∀𝑖, 𝑗, де 𝜀 — задане мале число. У даному випадку було вибрано 𝜀 ≈ 10−8. На Рис. 6Рис. 6 показанi початкова конфiгурацiя дна та розрахована в гiдродинамiчнiй задачi форма вiльної поверхнi, яку стацiонарний потiк утримує над нерозмивним па- горбом. З Рис. 6Рис. 6б випливає, що над вигнутим дном на поверхнi води утворюється вирва. Але, оскiльки в даному разi висота пагорба складає лише 10% вiд глибини води, ви- кривлення вiльної поверхнi виявляється дуже слабким — порядку 2 мм. Здатнiсть за- стосованої тут центральної KNP-схеми описувати такi незначнi змiни параметрiв течiї 144 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. а б Рис. 6. Початков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ни характеристик водного потоку. На Рис. 7Рис. 7 i 88 представлен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ях [1313,1414,2424]. Вiдсутнiсть нефiзичних осциляцiй у розв’язках свiдчить про те, що розроблена чисельна схема забезпечує акуратне моделювання ево- люцiї донної поверхнi як при слабкiй, так i при iнтенсивнiй взаємодiї наносiв з гiдроди- намiчною течiєю. Для оцiнки стiйкостi чисельного алгоритму було проведено довготривале моделю- вання розмиву пагорба в широкому каналi (розмiри розрахункової областi 2 × 2 км). Рис. 9Рис. 9 iлюструє отриману в цих розрахунках форму дна для 𝐴𝑔 = 0.001 при 𝑡 = 600 го- дин, а також вiдносну похибку щодо збереження об’єму наносiв, яка оцiнювалася за формулою (17)(17). Незважаючи на те, що з початку процесу минув достатньо великий час, результуюча похибка розрахунку не перевищила 2%. Геометричним параметром, що характеризує еволюцiю пiщаного пагорба у пото- цi води, є кут розмиву, утворений центральною вiссю нерiвностi й дотичною до неї в 145 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. а б в Рис. 7. Еволюцiя пагорба при 𝐴𝑔 = 0.001 (розмив нестацiонарною течiєю): а — 𝑡 = 25 год, б — 𝑡 = 50 год, в — 𝑡 = 100 год а б в Рис. 8. Еволюцiя пагорба при 𝐴𝑔 = 1 (розмив нестацiонарною течiєю): а — 𝑡 = 50 c, б — 𝑡 = 300 c, в — 𝑡 = 600 c 146 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. Рис. 9. Форма пагорба при 𝑡 = 600 год i вiдносна похибка морфологiчного розрахунку Рис. 10. Кут розмиву пiщаного пагорба при 𝐴𝑔 = 0.001 крайнiй точцi сформованої на днi зiркоподiбної структури. У роботi [2424] отримано наб- лижену аналiтичну оцiнку для кута розмиву 𝛼 пiщаного пагорба при слабкiй взаємодiї придонного потоку наносiв з течiєю води (𝐴𝑔 < 0.01): tg𝛼 = 3 √︀ 3(𝑚 + 1) 9𝑚− 1 . (18) При 𝑚 = 3, як було прийнято в наших розрахунках, iз формули (18)(18) одержуємо 𝛼 = arctg (3 √ 3/3) ≈ 22∘. На Рис. 10Рис. 10 представленi розрахованi контури донної поверхнi на рiвнi 1.5 см у рiз- нi моменти часу. Вони iлюструють динамiку розмиву пiщаного пагорба при слабкiй взаємодiї наносiв i течiї води. Як видно з графiкiв, еволюцiя пагорба вiдбувається в рамках трикутника, кут при вершинi якого (кут розмиву) близький до аналiтичного прогнозу [2424] i чисельних даних, представлених у статтях [1313,1414]. В ц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 кар’єри iстотно впливають як на гiдроди- намiчнi характеристики течiї, так i на морфологiю русла, що може призводити до пош- кодження та руйнування рiчкових iнженерних споруд — мостiв, пiдводних комунiкацiй, дамб обвалування та iн. 147 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. Область моделювання має вигляд прямокутного каналу, на днi якого розташований пiщаний кар’єр. Вважаємо, що момент 𝑡 = 0 вiдповiдає закiнченню розробки кар’єру, пiсля чого вiн починає iснувати за законами гiдроморфодинамiки. Схема розрахун- кової областi у вертикальнiй площинi 𝑂𝑥𝑧 i початковi умови задачi представленi на Рис. 11Рис. 11. Зазначимо, що для зручностi iнтерпретацiї результатiв у цьому випадку вiсь 𝑂𝑥 пов’язано з незбуреним дном. Характернi геометричнi параметри заглиблення — дiаметр 𝐷 i максимальна глибина 𝑑0. Його положення в площинi 𝑂𝑥𝑦 визначається координатою центральної точки (𝑋0, 𝑌0), так що на початку процесу отвiр кар’єру ло- калiзований в областi [𝑋1, 𝑋2]× [𝑌1, 𝑌2], де 𝑋1 = 𝑋0−𝐷/2; 𝑋2 = 𝑋0+𝐷/2; 𝑌1 = 𝑌0−𝐷/2; 𝑌2 = 𝑌0 + 𝐷/2. Будемо вважати, що заглиблення має форму конуса. Тодi при 𝑡 = 0 топографiю дна можна описати функцiєю 𝑧𝑏(𝑥, 𝑦, 0) = ⎧⎪⎨⎪⎩ −𝑑0 sin2 (︂ 𝜋(𝑥−𝑋1) 𝐷 )︂ sin2 (︂ 𝜋(𝑦 − 𝑌1) 𝐷 )︂ при (𝑥, 𝑦) ∈ [𝑋1, 𝑋2] × [𝑌1, 𝑌2]; 0 в iнших випадках. Якщо 𝑤0 — рiвень незбуреної вiльної по- Рис. 11. Конфiгурацiя розрахункової областi в центральному перерiзi площини 𝑂𝑥𝑧 верхнi, а 𝑄 [м3/с] – витрата води в каналi, то початковi умови для гiдродинамiчних параметрiв течiї мають наступнiй вигляд: ℎ(𝑥, 𝑦, 0) = 𝑤0 − 𝑧𝑏(𝑥, 𝑦, 0), 𝑢(𝑥, 𝑦, 0) = 𝑄 ℎ(𝑥, 𝑦, 0) , 𝑣(𝑥, 𝑦, 0) = 0. Як i в попереднiй задачi про розмив пiщаного пагорба, у постановку граничної задачi слiд включити вимоги фiксованостi рiвня й витрати води у вхiдному перерiзi каналу, нульових градiєнтiв усiх величин на виходi, а також умову вiдбиття потокiв води й наносiв на бiчних стiнках каналу. Розрахункова область, яка має розмiри 4000 × 2000 м, дискретизується вiдрiзками ∆𝑥 = ∆𝑦 = 20 м. Крок за часом обчислюється за критерiєм (6)(6) при 𝐶𝐹𝐿 = 0.5. При мо- делюваннi покладається властива для рiвнинних рiчок слабка взаємодiя мiж потоками води й наносiв (𝐴𝑔 = 0.001 i 𝑚 = 3 у формулi (5)(5)) при пористостi ґрунту 𝑝𝑚 = 0.4. При розв’язаннi гiдродинамiчних рiвнянь нахил тертя враховується з коефiцiєнтом Маннiнга 𝑛 = 0.01с/м1/3. Початковi умови для повної гiдроморфологiчної проблеми, в рамках якої розрахову- ються еволюцiя донної поверхнi й зумовленi нею змiни гiдродинамiчної течiї, визнача- ються з розв’язання гiдродинамiчних рiвнянь при фiксованому днi до набуття потоком стану рiвноваги. На Рис. 12Рис. 12 представленi результати моделювання рiвноважної течiї в каналi з початковим рiвнем води 𝑤0 = 5 м i витратою води 𝑄 = 5м3/с/ На днi каналу розташований пiщаний кар’єр дiаметром 𝐷 = 200 м, глибиною 𝑑0 = 5 м, центр якого 148 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. а б в г Рис. 12. Форма дна i рiвноважнi характеристики потоку над западиною: а — iзолiнiї донної поверхнi 𝑧𝑏, б — iзолiнiї поверхнi води 𝑤, в — розподiл поздовжньої швидкостi 𝑢, г — розподiл поперечної швидкостi 𝑣 має координати 𝑋0 = 400 м, 𝑌0 = 1000 м (див. Рис. 12Рис. 12а). На Рис. 12Рис. 12б, який iлюструє отриману в розрахунках форму вiльної поверхнi, можна бачити невеликий пiдйом рiвня води над заглибленням вiдносно початкового значення. Рис. 12Рис. 12в 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 його глибиною ℎ0 i швидкiстю над незбуреним дном 𝑈 (див. Рис. 11Рис. 11). На- приклад, представленi на Рис. 12Рис. 12 результати розрахункiв отримано при ℎ0 = 𝑤0 = 5 м, 𝑈 = 𝑄/ℎ0 = 1 м/с. Еволюцiю кар’єру в потоцi з цими початковими умовами iлюструє Рис. 13Рис. 13, на якому наведенi iзолiнiї функцiї 𝑧𝑏(𝑥, 𝑦) у рiзнi моменти часу. З цих графiкiв видно, що пiд впливом потоку води локалiзоване заглиблення трансформується в склад- ну систему ям i невеликих насипiв, яка дрейфує вниз за течiєю й поширюється вздовж поперечного перерiзу каналу. З якiсної точки зору картина розмиву пiщаного кар’єру виглядає наступним чином. Доннi наноси, якi транспортуються водною течiєю вище вiд 149 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. а б в г Рис. 13. Еволюцiя пiщаного кар’єру на днi каналу (розподiл функцiї 𝑧𝑏) при ℎ0 = 5 м, 𝑈 = 1 м/с, 𝐷 = 200 м, 𝑑0 = 5 м: а — 𝑡 = 100 год, б — 𝑡 = 300 год, в — 𝑡 = 500 год, г — 𝑡 = 700 год а б Рис. 14. Змiна форми вiльної поверхнi в процесi еволюцiї кар’єру (розподiл функцiї 𝑤) при ℎ0 = 5 м, 𝑈 = 1 м/с, 𝐷 = 200 м, 𝑑0 = 5 м: а — 𝑡 = 100 год, б — 𝑡 = 500 год 150 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. кар’єру, пiдходячи до нього, скочуються по верховому укосовi. Таким чином, iз цього боку кар’єр поступово засипається. Водночас, удар потоку води об низовий укiс викли- кає його розмив, через що центральна частина кар’єру сильно видовжується. Утворення бiчних ям пов’язане з появою над заглибленням поперечної складової швидкостi водної течiї, як показано на Рис. 12Рис. 12г. Ц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 на Рис. 14Рис. 14 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ю 𝑑max(𝑡), значення якої вiдповiдає максимальнiй глибинi кар’єру в заданий момент часу, та вели- чину 𝑥(𝑑max) — поздовжню координату точки, в якiй досягається максимальна глибина. Очевидно, що за функцiєю 𝑥(𝑑max(𝑡)) можна охарактеризувати швидкiсть, з якою за- глиблення перемiщується вздовж русла. Рис. 15Рис. 15 iлюструє вплив геометричних i гiдродинамiчних параметрiв розглянутої за- дачi на еволюцiйнi характеристики кар’єру. Представленi на Рис. 15Рис. 15а залежностi 𝑑max(𝑡) i 𝑥(𝑑max) отримано для рiзних початкових глибин кар’єру 𝑑0 при фiксованих його дiа- метрi й характеристиках водної течiї — 𝐷 = 200 м, ℎ0 = 5 м, 𝑈 = 1 м/с. З цих графiкiв випливає, що градiєнт функцiї 𝑑max(𝑡) для глибоких ям набагато вищий, нiж для мiлких, тобто, глибокi кар’єри засипаються швидше вiдносно своєї початкової глибини. З iншо- го боку, судячи з характеру отриманих залежностей 𝑥(𝑑max), вони значно повiльнiше пересуваються по дну. Вплив швидкостi води 𝑈 на моделювання еволюцiї розмиву кар’єрiв зумовлений використанням степеневої формули (5)(5) для опису потоку наносiв. Зокрема, при 𝑚 = 3 отримуємо кубiчну залежнiсть потужностi наносiв вiд швидкостi води, з якої випливає, що при 𝑈 < 1доннi процеси iстотно сповiльнюються, а при 𝑈 > 1 — прискорюються. Цей висновок пiдтверджується Рис. 15Рис. 15б, на якому зображенi графiки функцiй 𝑑max(𝑡) i 𝑥(𝑑max), отриманi при фiксованих характеристиках процесу 𝑑0 = 12 м, 𝐷 = 200 м, ℎ0 = 5 м для рiзних значень швидкостi води в каналi — 𝑈 = 0.5, 1 i 2 м/с Залежнiсть морфологiчних процесiв вiд глибини води в каналi ℎ0 не носить такого явного характеру, як вiд швидкостi 𝑈 . З представлених на Рис. 15Рис. 15в графiкiв випли- ває, що пiщанi кар’єри в глибоких акваторiях засипаються набагато довше i дрейфу- ють значно повiльнiше, нiж у мiлких. Цi результати було отримано при 𝑑0 = 10 м, 151 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. а б в Рис. 15. Еволюцiйнi характеристики кар’єру 𝑑max(𝑡) i 𝑥(𝑑max): а — вплив початкової глибини кар’єру, б — вплив швидкостi води в каналi, в — вплив глибини води в каналi 𝐷 = 200 м, 𝑈 = 1 м/с. Вiдзначимо також практично лiнiйну залежнiсть еволюцiйних характеристик кар’єру вiд глибини води в каналi. Розрахунки показують, що еволюцiя кар’єру вiдбувається в рамках трикутника, по- будованого дотичними до контурiв функцiї 𝑧𝑏(𝑥, 𝑦) одного й того малого ж рiвня (див. Рис. 16Рис. 16). Кут 𝛼 при вершинi цього трикутника (кут розмиву) є одним з головних геомет- 152 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. ричних параметрiв, якi характеризують еволюцiю пiщаного кар’єру в рiчковому потоцi, оскiльки вiн дозволяє оцiнити поширення форми дна в обох напрямках — поперечному i поздовжньому. На Рис. 17Рис. 17 представленi розрахованi залежностi кута розмиву 𝛼 вiд початкової глибини кар’єру 𝑑0 i швидкостi води 𝑈 . Тут можна бачити, що кут розмиву зростає зi збiльшенням глибини кар’єру нелiнiйно. Ймовiрно, це зумовлено наявнiстю багатьох факторiв, якi впливають на процес. Залежнiсть 𝛼 вiд швидкостi потоку 𝑈 Рис. 16. Кут розмиву 𝛼 пiщаного кар’єру при 𝑑0 = 5 м, 𝐷 = 200 м, 𝑈 = 1 м/с, ℎ0 = 5 м носить неоднозначний характер. Найбiль- шi кути розмиву в розглянутому дiапазонi швидкостей вiдповiдають 𝑈 = 0.5 м/с. Цей результат пiдтверджує, що при 𝑈 < 1 до- мiнує поширення форми дна в поперечно- му напрямку, яке завершується утворенням великих ерозiйних зон упоперек потоку. В цiлому з графiкiв, наведених на Рис. 17Рис. 17, випливає, що при реалiстичних параметрах течiї кути розмиву великих донних заглиб- лень лежать у межах вiд 22∘ до 34∘. На Рис. 18Рис. 18 представлено вiдносну по- хибку збереження об’єму наносiв (17)(17), от- риману при моделюваннi розмиву кар’єрiв рiзної початкової глибини. Як видно з гра- фiка, у всiх випадках вона не перевищувала 1%. Цей результат разом з вiдсутнiстю нефiзичних осциляцiй при оцiнюваннi форм дон- ної поверхнi (див. Рис. 13Рис. 13 i 1616) свiдчить про те, що запропонований алгоритм забезпечує акуратне моделювання транспорту придонних наносiв. Окрiм того, завдяки застосуван- ню зв’язаного пiдходу при розв’язаннi системи рiвнянь (1)(1), (2)(2), вiн здатен забезпечити розрахунок оберненого впливу морфологiчних процесiв на течiю води. Рис. 17. Залежностi кута розмиву 𝛼 пiщаного кар’єру вiд початкової глибини 𝑑0 при рiзних швидкостях води в каналi (𝐷 = 200 м, ℎ0 = 5 м) Рис. 18. Вiдносна похибка розрахункiв для кар’єрiв рiзної глибини (𝑈 = 1 м/с, ℎ0 = 5 м) 153 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. 6. ВИСНОВКИ 1. Розроблено математичну модель i чисельний алгоритм для розрахунку гiдромор- фодинамiчних процесiв у вiдкритих водоймах, який поєднує центральну схему другого порядку Курганова—Ноелля—Петрової (KNP-scheme) для iнтегрування гiдродинамiчних рiвнянь зi схемою «проти течiї» п’ятого порядку, за якою роз- раховується еволюцiя донної поверхнi. Остання ґрунтується на зваженiй iстотно неосцилюючiй апроксимацiї (WENO-approximation) потоку наносiв. 2. Тестування алгоритму, проведене на в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 часов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д 20∘ до 35∘ i залежать вiд їхньої початкової глибини та швидкостi руслового потоку. 5. Зробленi оцiнки важливi для прогнозування й попередження негативних техно- генних i природних явищ в рiчкових акваторiях, якi можуть бути викликанi ди- намiкою великих донних утворень та їхньою взаємодiєю з гiдротехнiчними спору- дами, берегами рiчок тощо. ЛИТЕРАТУРА [1] Wu W. Computational river dynamics. –– London : Taylor and Francis, 2008. –– 487 p. [2] Castro-Diaz M. J., Fernandez-Nieto E. D., Ferreiro A. M. Sediment transport models in shallow water equations and numerical approach by high order finite volume methods // Computers & Fluids. –– 2008. –– Vol. 37. –– P. 299–316. [3] Горбань I. М. Чисельне моделювання еволюцiї нерiвностей великого масштабу на рiчковому днi // Прикладна гiдромеханiка. — 2015. — Т. 17(89), № 1. — С. 21–36. [4] Hudson J., Sweby P. K. Formulations for numerically approximating hyperbolic systems governing sediment transport // Journal of Scientific Computing. –– 2003. –– Vol. 19. –– P. 225–252. 154 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. [5] One-dimensional model for transient flows involving bed-load sediment transport and changes in flow regimes / L. Goutiere, S. Soares-Frazao, C. Savary et al. // Journal of Hydraulic Engineering. –– 2008. –– Vol. 134, no. 6. –– P. 726–735. [6] Kubatko E., Westerink J., Dawson C. An unstructured grid morphodynamic model with a discontinuous galerkin method for bed evolution // Ocean Modelling. –– 2006. –– Vol. 15. –– P. 71–89. [7] Long W., Kirby J. T., Shao Z. A numerical scheme for morphological bed level calcula- tions // Coastal Engineering. –– 2008. –– Vol. 55. –– P. 167–180. [8] Сухомел Г. И. Исследование гидравлики открытых русел и сооружений. — Киев : Наукова думка, 1965. — 112 с. [9] Liu X.-D., Osher S., Chan T. Weighted essentially non-oscillatory schemes // Journal of Computational Physics. –– 1994. –– Vol. 115. –– P. 200–211. [10] Kurganov A., Noelle S., Petrova S. Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton—Jacobi equations // SIAM Journal on Scientific Com- puting. –– 2001. –– Vol. 23, no. 3. –– P. 707–740. [11] Kurganov A., Lin C.-T. On the reduction of numerical dissipation in central-upwind schemes // Communications in Computational Physics. –– 2007. –– Vol. 2, no. 1. –– P. 141–163. [12] Gottlieb S., Shu C. W., Tadmor E. Strong stability-preserving high order time dis- cretization methods // SIAM Review. –– 2001. –– Vol. 43. –– P. 89–112. [13] Hudson J., Sweby P. K. A high-resolution scheme for the equations governing 2d bed- load sediment transport // International Journal for Numerical Methods in Fluids. –– 2005. –– Vol. 47. –– P. 1085–1091. [14] Benkhaldoun F., Sahmim S., Seäıd M. A two-dimensional finite volume morphodynamic model on unstructured triangular grids // International Journal for Numerical Methods in Fluids. –– 2010. –– Vol. 63. –– P. 1296–1327. [15] Kurganov A., Tadmor E. New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations // Journal of Computational Physics. –– 2000. –– Vol. 160. –– P. 241–282. [16] Горбань В. О., Горбань I. М. Чисельнi моделi нестацiонарних руслових процесiв // Прикладна гiдромеханiка. — 2013. — Т. 15(87), № 4. — С. 19–39. [17] Yee H. C. Construction of explicit and implicit symmetric TVD schemes and their applications // Journal of Computational Physics. –– 1987. –– Vol. 68. –– P. 151–179. [18] Kurganov A., Levy D. Central-upwind schemes for the Saint-Venant system // ESAIM: Mathematical Modelling and Numerical Analysis. –– 2002. –– Vol. 36. –– P. 397–425. 155 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. [19] Johnson H. K., Zyserman J. A. Controlling spatial oscillations in bed level update schemes // Coastal Engineering. –– 2002. –– Vol. 46. –– P. 109–126. [20] Jensen J. H., Madsen E. Ø., Fredsøe J. Oblique flow over dredged channels: II. Sediment transport and morphology // Journal of Hydraulic Engineering. –– 1999. –– Vol. 125, no. 11. –– P. 1190–1198. [21] Hurten A. High resolution schemes for hyperbolic conservation laws // Journal of Com- putational Physics. –– 1983. –– Vol. 49. –– P. 357–365. [22] LeVeque R. J. Balancing source terms and flux gradients in high-resolution Godunov methods: The quasi-steady wave-propagation algorithm // Journal of Computational Physics. –– 1998. –– Vol. 146. –– P. 346–356. [23] Hing Y., Shu C.-W. High order finite difference WENO schemes with the exact conser- vation property for the shallow water equations // Journal of Computational Physics. –– 2005. –– Vol. 208. –– P. 206–227. [24] Hudson J. Numerical techniques for morphodynamic modelling : Ph. D. thesis / J. Hud- son ; University of Reading. –– 2001. –– 185 p. REFERENCES [1] W. Wu, Computational river dynamics. London: Taylor and Francis, 2008. [2] M. J. Castro-Diaz, E. D. Fernandez-Nieto, and A. M. Ferreiro, “Sediment transport models in shallow water equations and numerical approach by high order finite volume methods,” Computers & Fluids, vol. 37, pp. 299–316, 2008. [3] I. M. Gorban, “Numerical simulation of evolution of large-scale irregularities on a river bottom,” Applied Hydromechanics, vol. 17(89), no. 1, pp. 21–36, 2015. [4] J. Hudson and P. K. Sweby, “Formulations for numerically approximating hyperbolic systems governing sediment transport,” Journal of Scientific Computing, vol. 19, pp. 225–252, 2003. [5] L. Goutiere, S. Soares-Frazao, C. Savary, T. Laraichi, and Y. Zech, “One-dimensional model for transient flows involving bed-load sediment transport and changes in flow regimes,” Journal of Hydraulic Engineering, vol. 134, no. 6, pp. 726–735, 2008. [6] E. Kubatko, J. Westerink, and C. Dawson, “An unstructured grid morphodynamic model with a discontinuous Galerkin method for bed evolution,” Ocean Modelling, vol. 15, pp. 71–89, 2006. [7] W. Long, J. T. Kirby, and Z. Shao, “A numerical scheme for morphological bed level calculations,” Coastal Engineering, vol. 55, pp. 167–180, 2008. [8] G. I. Suhomel, A study of hydraulics of open channels and structures. Kiev: Naukova Dumka, 1965. 156 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. [9] X.-D. Liu, S. Osher, and T. Chan, “Weighted essentially non-oscillatory schemes,” Journal of Computational Physics, vol. 115, pp. 200–211, 1994. [10] A. Kurganov, S. Noelle, and S. Petrova, “Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton—Jacobi equations,” SIAM Journal on Scientific Computing, vol. 23, no. 3, pp. 707–740, 2001. [11] A. Kurganov and C.-T. Lin, “On the reduction of numerical dissipation in central-upwind schemes,” Communications in Computational Physics, vol. 2, no. 1, pp. 141–163, 2007. [12] S. Gottlieb, C. W. Shu, and E. Tadmor, “Strong stability-preserving high order time discretization methods,” SIAM Review, vol. 43, pp. 89–112, 2001. [13] J. Hudson and P. K. Sweby, “A high-resolution scheme for the equations governing 2d bed-load sediment transport,” International Journal for Numerical Methods in Fluids, vol. 47, pp. 1085–1091, 2005. [14] F. Benkhaldoun, S. Sahmim, and M. Seäıd, “A two-dimensional finite volume morphodynamic model on unstructured triangular grids,” International Journal for Numerical Methods in Fluids, vol. 63, pp. 1296–1327, 2010. [15] A. Kurganov and E. Tadmor, “New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations,” Journal of Computational Physics, vol. 160, pp. 241–282, 2000. [16] V. O. Gorban and I. M. Gorban, “Numerical models of non-stationary river processes,” Applied Hydromechanics, vol. 15(87), no. 4, pp. 19–39, 2013. [17] H. C. Yee, “Construction of explicit and implicit symmetric TVD schemes and their applications,” Journal of Computational Physics, vol. 68, pp. 151–179, 1987. [18] A. Kurganov and D. Levy, “Central-upwind schemes for the Saint-Venant system,” ESAIM: Mathematical Modelling and Numerical Analysis, vol. 36, pp. 397–425, 2002. [19] H. K. Johnson and J. A. Zyserman, “Controlling spatial oscillations in bed level update schemes,” Coastal Engineering, vol. 46, pp. 109–126, 2002. [20] J. H. Jensen, E. Ø. Madsen, and J. Fredsøe, “Oblique flow over dredged channels: II. Sediment transport and morphology,” Journal of Hydraulic Engineering, vol. 125, no. 11, pp. 1190–1198, 1999. [21] A. Hurten, “High resolution schemes for hyperbolic conservation laws,” Journal of Computational Physics, vol. 49, pp. 357–365, 1983. [22] R. J. LeVeque, “Balancing source terms and flux gradients in high-resolution Godunov methods: The quasi-steady wave-propagation algorithm,” Journal of Computational Physics, vol. 146, pp. 346–356, 1998. 157 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. [23] Y. Hing and C.-W. Shu, “High order finite difference WENO schemes with the exact conservation property for the shallow water equations,” Journal of Computational Physics, vol. 208, pp. 206–227, 2005. [24] J. Hudson, Numerical techniques for morphodynamic modelling. phdthesis, University of Reading, 2001. И. Н. Горбань Моделирование динамики песчаных карьеров в русловых потоках Развит численный метод для расчета совместной эволюции водного течения и дон- ной поверхности в открытых водоемах, учитывающий взаимное влияние гидро- динамического и морфологического процессов. Он сочетает центральную схему «против течения» Курганова—Ноелля—Петровой для интегрирования уравнений мелкой воды и взвешенную существенно неосциллирующую (WENO) схему пято- го порядка, по которой рассчитывается эволюция донной поверхности. Мощность придонного потока наносов оценивается по степенной формуле Грасса. Модели- рование тестовой задачи о размыве потоком воды конического песчаного холма выявило долговременную устойчивость полученного решения, а также его соот- ветствие данным других авторов. Метод применен к расчету эволюции глубоких песчаных карьеров на речном дне. Полученные результаты демонстрируют форми- рование на донной поверхности сложного рельефа, когда локализованное углубле- ние трансформируется в систему ям и небольших насыпей, перемещающихся вслед за течением. В то же время, углубление расширяется поперек потока, что приво- дит к формированию протяженных эрозионных зон на речном дне. Установлено, что величина угла размыва, характеризующего эволюцию карьера, зависит от на- чальной глубины карьера и параметров потока воды и находится в пределах от 20∘ до 35∘. Полученные оценки «живучести» песчаных карьеров указывают на нелинейную зависимость интенсивности их засыпания и скорости перемещения от характеристик течения. Моделирование эволюции свободной поверхности выяви- ло, что над неровностью дна генерируется интенсивный волновой фронт, форма которого изменяется в соответствии с процессами, происходящими на донной по- верхности. Найденные оценки могут быть использованы при организации меропри- ятий по предотвращению негативных техногенных и природных явлений в речных акваториях, вызванных взаимодействием крупномасштабных донных образований с гидротехническими сооружениями и берегами. КЛЮЧЕВЫЕ СЛОВА: морфологический процесс, донная поверхность, транспорт наносов, песчаный карьер, угол размыва I. M. Gorban Modeling of sandpit dynamics in river flows A two-dimensional numerical method for simulating of coupled evolution of the water flow and the bed in open reservoirs is developed considering the mutual influence of the morphological processes. It combines the Kurganov-Noelle-Petrova central upwind scheme for integration of shallow water equations and the fifth-order Weighted Essen- tially Non-Oscillatory (WENO) Scheme to describe the evolution of bed surface. The 158 ISSN 2616-6135. ГIДРОДИНАМIКА I АКУСТИКА. 2018. Том 1(91), № 2. С. 132132–159159. bedload sediment transport rate is estimated by power-law Grass formula. Simula- tion of the test problem on erosion of a conical sand dune in a water flow revealed the longterm stability of the obtained solution, as well as its relevance to the data presented by other researchers. The method was used to calculate the evolution of deep sand quarries on the river bed. The obtained results demonstrate formation of a compound relief form on the bed when the localized hollow is transformed into the system of pits and small mounds moving with the flow. At the same time, the deep- ening expands across the flow that leads to the formation of extended erosion zones on the river bottom. The angle of spread characterizing the evolution of the quarry is found to depend on the initial quarry depth and water flow parameters and ranges from 20∘ to 35∘. The obtained estimates of “survivability” of sand quarries indicate the nonlinear dependence of both the intensity of their filling and motion velocity on the flow parameters. Simulation of free surface evolution revealed an intense wave front generation over the unevenness of the bottom. The shape of this front changes in ac- cordance with the processes occurring at the riverbed. The obtained estimates may be used when organizing the activities to prevent the negative technological and natural phenomena in river water areas initiated by interaction of large-scale bed roughness with hydraulic structures and banks. KEY WORDS: morphological process, bed surface, sediment transport, sand quarry, angle of spread 159 ВСТУП ОСНОВНІ РІВНЯННЯ ЧИСЕЛЬНА СХЕМА Моделювання гідродинамічної течії Дискретизація морфологічного рівняння ТЕСТОВІ РОЗРАХУНКИ Малі збурення вільної поверхні над неплоским нерозмивним дном Розмив піщаного пагорба в стаціонарному потоці Розмив піщаного пагорба в нестаціонарному потоці РОЗМИВ ПІЩАНИХ КАР'ЄРІВ ВИСНОВКИ

Моделювання динаміки піщаних кар'єрів у руслових потоках

Репозитарії

Схожі ресурси