Математическое моделирование волнового движения воды в узко-глубоком непризматическом водохранилище с учётом перелива (перехлёста) воды через гребень плотины

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

Том 19, № 2, 2014

Математическое моделирование волнового движения

воды в узко-глубоком непризматическом водохранилище с учётом перелива (перехлёста) воды

через гребень плотины

И. Д. МУЗАЕВ1'2, Н.И. МУЗАЕВ2 1 Северо-Осетинский государственный университет им. К.Л. Хетагурова, 2Центр геофизических исследований ВНЦ РАН и РСО-А, Владикавказ, Россия

e-mail: muzaevid@mail.ru

Музаев И.Д., Музаев Н.И. Математическое моделирование волнового движения воды в узко-глубоком непризматическом водохранилище с учётом перелива (перехлёста) воды через гребень плотины // Вычисл. технологии. 2014. Т. 19, № 2. С. 94-106.

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

Ключевые слова: водохранилище, гребень плотины, перелив, водохранилище уз-коканьонного типа.

Muzaev I.D., Muzaev N.I. Mathematical modeling of wave motion of water in a narrow-deep non prismatic reservoir with the water overflow (overlap) over the dam // Comput. Technologies. 2014. Vol. 19, No. 2. P. 94-106.

The boundary value problem of water wave motion in a narrow-canyon type reservoir, when the waves are formed after landslide or avalanche invasion has been formulated and solved. Overflow over the dam was accounted for the first time in the history of scientific research on this problem.

Key words: tailings reservoir, the dam, overflow, narrow-canyon type reservoir.

Наличие глубоких и обширных водохранилищ в горных и предгорных местностях ставит перед исследователями ряд актуальных задач, связанных с образованием высоких волн в случае вторжения обвально-оползневого массива горной породы в чашу водохранилища. Правильная постановка и решение таких задач позволяют прогнозировать и предотвращать либо смягчать те последствия и ущерб, которые может вызвать образование разрушительных волн. Имеется много публикаций, касающихся постановки и решения конкретных задач, моделирующих вышеуказанный процесс. При этом предлагаются решения как аналитическими, так и конечно-разностными и конечно-элементными методами. В работах [1-3] задачи поставлены и решены в рамках линейной теории поверхностных волн малой амплитуды. Чаша водохранилища схематизирована в виде ограниченного или полуограниченного прямоугольного параллелепипеда.

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

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

Предположим, что в прямоугольной системе координат хОуг часть пространства, ограниченная условиями

0 < х < Ь, -

В(х, г) 2

< У <

В(х, г) ~2 :

-Н<г< 0,

представляет водохранилище узкоканьонного типа. Ось Ох направлена вдоль каньона, Оу — поперёк, Ог — вертикально вверх. В створе, ограниченном условиями

х

0,

-В

< У <

2

2

-Н<г< 0,

помещена плотина. Здесь Ь — длина водохранилища, В(х, г) — ширина каньона, зависящая от продольной координаты х и вертикальной координаты г, Н — глубина воды.

Рассмотрим волновое движение воды, вызванное вторжением массива горной породы с некоторого участка борта в данные заполненное водохранилище.

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

д2р д2р дх2 дг2

1 дВ(х, г) др В (х,г) дх дх В(х,г) дг

1 дВ

др дх

др

х=0

х=Ь

д2р др

дг2

дг

0,

др дг

х=0

г=-Н

др

Р = д* = 0 при г

0 и г = 0,

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

0

0

Функция д(х,г,£) берется в виде

д

I 1 при х0 — а < х < хо + а, 5дх) = \

I 0 при 0 ^ х ^ х0 — а и х0 + а < х < Ь,

(м при 0 <г<го, , .

ф) = \ п _ , (7) I 0 при I > 10,

где х0 — абсцисса центра фронта, 2а — ширина фронта обвально-оползневого массива, и0 — скорость вторжения, ¿0 — продолжительность времени вторжения. Граничное условие на напорной грани плотины

Эр dx

= v

x=0

учитывает перелив (перехлёст) воды через гребень плотины (рис. 1).

Легко заметить, что непризматическая конфигурация водохранилища смоделирована изменением ширины каньона B(x,z). Она зависит как от продольной координаты x, так и от вертикальной координаты z.

В связи с тем что коэффициенты дифференциального уравнения (1) зависят от пространственных координат, решение начально-краевой задачи (1)-(7) аналитическими методами связано с большими математическими трудностями. Эти трудности легко преодолеваются, если непризматическую конфигурацию водохранилища аппроксимировать экспоненциальными функциями следующего вида:

B (x,z) = BoeSlX eS2Z,

Bo = const, Si > 0, S2 > 0,

где B0 — длина гребня плотины, параметр S1 характеризует интенсивность расширения водохранилища в плане, параметр S2 характеризует интенсивность сужения водохранилища (каньона) по глубине.

Рис. 1. Схема процесса перелива (перехлёста) воды через гребень плотины

Подставив выражения (5), (6), (7) и (9) в (1)-(4), получим

д2р д2р др др дх2 дг2 дх дг

р

др

Ж при*

д1(х)и(г)^2(г) ^х

Во 6 6 ,

0, г = 0,

др дх

.ММ;), |

х=о

х=Ь

д2р др

дг2

д

дг

0,

z=0

др дг

0,

z=—Н

где .(г) — скорость перелива (перехлёста) воды через гребень плотины, (г) и ¿2(г) — вспомогательные функции, обеспечивающие гладкое сопряжение скоростей вторжения обвально-оползневого массива и скорости перелива с их нулевыми значениями. Скорость перелива воды через гребень плотины .(г) зависит от возвышения уровня свободной поверхности воды у плотины над её гребнем:

.

9

при — е < г < 0, при — Н ^ г ^ —е,

?(г) = Мп)л/2дп3/2(г)Во,

здесь д(г) — расход перелива воды через гребень плотины, ^(п) — коэффициент расхода, Во — длина плотины над отметкой её гребня. Выражение (15) представляет известную формулу водослива. Коэффициент расхода ^ зависит от величины п(*):

Мп)

0.402 при п(г) > 0, 0 при п(г) ^ 0.

Выражения (14) и (15) показывают, что расход перелива воды через гребень плотины распределён (рассредоточен) на участке плотины — е < г < 0, где е — бесконечно малая величина, т. е. е ^ 0. С точки зрения быстрой сходимости тригонометрических рядов, которые получаются в процессе решения начально-краевой задачи, целесообразно сглаживающие функции ¿1 (г) и ¿2(г) представить в виде

¿1 (г) = г3 (г + е)3 А, ^2(г) = г 3(г, + ^)3 В,

е6 к6

где к — мощность обвально-оползневого массива, постоянные А и В определяются из условий

— £ — Н

140

е6 к6 .

¿1(г)йг = 1, ¿2(г)^г = 1, А = —, В

оо

В результате применения интегрального преобразования Лапласа по времени г и

подстановки

р

(Ь — х)3 3Ь2

0

е

0

относительно новой неизвестной функции *(х,г) начально-краевая задача (10)-(13) принимает следующий вид:

д2* + д2* + — д* -2Ф - ) ^ + + 52 дг — 1"* = /(х,г)е 2х,

где

д* 51 1

х=о

2 г д*

р * + д^х

д* 51 1

дх — I1 *

д* дг

х=Ь

z=0

z=—H

—-^

/(х,г) = -(Р) (Ь — х)2

¿1(г)2(Ь^ + ¿'/(г)(Ь — х)3

Ь2

+ ЗД

Ь2

(Ь — х)3' 3Ь2

3Ь2

1

+ ^^(^(р) — е-* хе—'^. Во

Далее вводятся новая функция Ф(х, г) и дифференциальный оператор О

Ф д* 51Ф п д2 д2 с д -2 Ф = —--1 *, О =--1---+ -2---

дх 2 дх2 дг2 2 дг 4

Приложим оператор О к функции ф

~ д ~ 5 ~ О[ - ] = ддхО[ Ф ] — -2 О[ Ф ].

В развёрнутом виде выражение (25) имеет вид д2ф д2 ф „ дф д / ^

дх2 дг2 2 дг 4

дх

+ ттт + ---^Ф = тт" е"Г"х/

-1

2

Граничные условия (21) и (22) относительно функции Ф записываются как

Ф

0, Ф

х=о

х=Ь

2Ф дФ рФ+ддг

z=0

О дФ

дг

0.

z=—Н

Решение дифференциального уравнения (26) при граничных условиях (27) находится в виде следующего тригонометрического ряда по синусам:

те

Ф(х, г) = ^ Фп(г) 81п а„.х, аг.

п=1

Лп^} ^п

ПП

Т'

Разложим правую часть уравнения (26) в ряд Фурье по синусам в интервале (0, Ь):

— (е 2 7(х, г) J — — е 2 7(х,г) = ^Фп(г) 81

—1 «1х

— е 2 х -

2

г) = ¿п(г) 81папх,

п=1

0

0

0

0

$п(г) = у

дх

вт ап хё,х

ь

^ У ^ап сое апх + у вт апх) е^/(х, г)^х. 0

В результате подстановки значения функции /(х, г) из (23) в (30) и вычисления интегралов для 5п(г) получается выражение

где

¿п(г) = — Ь (г)—

ап + 4"

1 — (—1)пе51 ь I —

2

3Ь

^(рЖ/(г)ап—

2апЬ

2 + 52

ап +4

+

25" ап

2 -ап + 4

г

1 — (—1)пе 51 ь

+

2

2 5

+ 2т2 ¿1 (г )«п + -2- ¿2(г)ВДе-^х

3Ь3

ЬВ0

ап = (Ь — х)3е 21 ( ап сое апх + —1 вт апх ) ^х,

5Т 2

Дп = I (х)е 2 х ( ап сое апх + вт апх ) ^х.

Подставляя выражения (29) и (30) в (26)-(28) и приравнивая коэффициенты при одинаковых синусах в правых и левых частях, получим

-^г^ + ^ — К + т)Ф п = ад,

^Ф п

Р2Ф п + д п

^г

^Ф п

х=0

^г

х=-Н

Решение дифференциального уравнения (34) с граничными условиями (35) имеет

вид

~ 1 [ ~

Фп

Ап .

где

р2 — д-21 в

С"

Ап сЬ АпН + -2 8Ь ЛпЯ |

ь

2

ь

ь

0

0

X

С2 = —-

д^пВ«

Лп еЬ ЛпН + 8Ь АпН)

7п

\

2 52

д I Лп — -4-1 8Ь ЛпН

-2 ' Ап еЬ ЛпН + — 8Ь ЛпН

—Н

В = — Лт/е 52'

52

Лп еЬ Лп

Фп(г)^г.

В результате подстановки значения Фп(г) из (31) в (41), вычисления интегралов и перехода к пределу при е ^ 0 получим следующее выражение для Вп:

Ф 1 52 Вп = — I Лп еЬ ЛпН + — 8Ь ЛпН

4

Ь2

ап

2 + ап +4

2—1 Т3

2Ьап 2—1ап

—с + , , о

-12

ап +

2 + ап +4

—/• \ I 2^пКп—/ \

Ф

ЛпЬВ0

где

- Н

Кп = I е-г2 ( Лп еЬ Лп(Н + г) + у8Ь Лп(Н + г) | ¿2(г)йг.

Подставив выражение (41) в (36), получим

Ф п(0,г) = —д

4ап

-2

Ь3 I ап + --

1 — ("1)"е*') Ж + дЛп1 х

Р2 + тп

Ь3

х

2Ьа п

-12

+

2-1 ап (1 — (—1)пе51 ^

ап +

-12

ап +

д(р) р2 + тГ

2ДпКп

Ф

, -2 \ Р2 + '

ЬВ0 I Лп еЬ ЛпН + — 8Ь ЛпН

2

Переходя от изображения к оригиналу, имеем

4ап

Фп(0,Ь) = —д

(1 — (—1)пе51 ь) 1 У д(т)в1п7п(* — т)^т+

52

ЬМ ап + ^

+дЛ

2-1

( \

2Ьап 2-тап

2 +

ап + 4

+

V

2 + 52

ап +4

— I д(т)в1п7п(^ — т)^т— 7п ] 0

2ДпКп

1

— и(т) в1п 7п(^ — т)^т. —2 \ 1п J

ЬВ0 ( Лп сЬ ЛпН + у вЬ ЛпН I 0

Далее обратным ходом сначала в результате решения дифференциального уравнения

р

п=1

где

/п

—1 . 2

— в1п апх + ап сов апх

е 2 х.

2 + —2 ап + -4-

Уравнение волновой поверхности определяется дифференцированием потенциала <^(х, 0,Ь) по переменной Ь [10-12]

п(х,Ь) = — 1 , п(х,Ь) = —1 с/(Ь) + 1 ^ Фп(Ь)1п(х).

п=1

Произвольная функция с/(Ь) находится из закона сохранения объёма жидкости в ь ь ь -Н

/ в°е-'?1Хч<х'г)А=—1ф )<гт+Л / ^^

0 0 0 0 0

Подставив выражение (47) в (48), имеем

С

д51

(е^1ь — 1)В0

д(т)^т —

д2аЛ/0(Ь) 5-

В0 е^ь - 1:

/0

Ь при 0 < Ь < Ь0, 1 при Ь > ¿0.

ь

2

ь

В итоге для закона колебания уровня воды у плотины получается выражение

п(г) = —

— 1

£

п=1 „2

-12

ап +

4 ап (1 — (—1)пе51 ^

Ь3

-12

ап + 4

4-1ап^2д - 2 4

ЬМ ап + ^

ь

а

п

4-2ап ^Тд -2

ЬМ ап + ^

ео87п(г — т) )п3/2(т+

2ак-1

В0

и(т)^т+

2ЯпКп / и(т) ео8 7п(г — т)^т

£

0

п=1 ЬВ^ Лп еЬ ЛпН + 8Ь ЛпН| ( аЩ + --

представляющее собой нелинейное интегральное уравнение типа уравнения Вольтерра относительно искомой функции п(0, ¿) = п(*), описывающей процесс колебания уровня воды у плотины. При -1 = -2 = 0

п(0,г) = —I — £ 4(1 ,2(—31)п) ео8 7п(г — т )\^(т )п3/2 (т +

Ь ^ апЬ3

п=1 п

2ак

+ТГ-Г I и(

В0Ь

У и(т)^т + £ ЬВ' и(т) ео8 7п(г — т)^т.

п п=1 п П

0 п=1 0 В результате вычисления интегралов выражение (50) можно привести к виду

ь

п

В0Ь В0 п= ап7п Ь ]

п—1 п

/0

г при 0 < г < г0, г0 при г > г0,

81п 7п*

при 0 < г < г0

/п(*)=<2 /, Ч . г0 .

2 ео8 7п (г — — I 81п 7п — при г > ¿0,

4

• /-1-тт ПП

ап = 7— ео8 апх0 81п апа, 7п = V дап ш апН, ап = —, Ьап Ь

2 ^ 2 1 - (-1)п К (г,т) = 1 + 7£вп ео8 7п(г — т),

п=1

Ь ап

ь

2

ь

ь

ь

Поскольку (50) является нелинейным интегральным уравнением, решение строгими аналитическими методами связано с большими математическими трудностями. Поэтому применяется численно аналитический метод решения. В качестве первого приближения использовано выражение

2аНп0 п0 ^ ап Ш апН

т = „ т т) + — у -/n(t),

Б0Ь Б0 ^ апуп

п= 1

второго

г

%(*) = т к (г,т Ыт )1ш(т )|3/2ж

и т. д.

о г

Цп(1) = Пп-г(г) \ к(г,тК-1 (т)|Пп—1 (т)|3/2^т.

Для вычисления интеграла

г

I к (г,т (т )|п1(т )|3/2^т о

использована формула трапеции г

I к (*,тМт)Ыт)Г ¿т « Агк 0)М1(0)|П1(0>|3/2 + к Е.Ф, МЫ«)|3/2.

г-1

+ £ к (1,3 а^Мз а^)Ыз А^)|3/2.

3 = 1

Скорость и расход перелива (перехлёста) воды через гребень плотины вычисляются формулами

г;(*) = /при ф) > 0.

4 ) \0 при ф) < 0,

= 1 |3/2 при п(*) > 0,

) |0 при ц(г) < 0.

В качестве примеров на рис. 2, 3 представлены графики зависимостей от времени уровня воды в створе плотины, объёма, расхода и скорости перелива (перехлёста) воды через гребень плотины при различных значениях входных параметров. Из рисунков следует, что с течением времени из-за перелива воды амплитуда образованных волн, расход и скорость перелива воды через гребень плотины уменьшаются. Из приведённых графиков также видно, что образованная в середине водохранилища волна добегает до плотины за ~ 100 с, проходя за это время расстояние 2500 м со скоростью 25 м/с.

Согласно классической теории волн для скорости волны имеется формула

^ <52)

С

Г], М

0 100 200 300 400 500 600 700 800 900 t, с 0 100 200 300 400 500 600 700 800 900 t, с

90 80 70 60 50 40 30 20 10 О

q, м3/с

V, м/с

-г-Т-Т-Т-т-

*  i

WJMJii

о 100 200 300 400 500 600 700 800 900 t, с 0 100 200 300 400 500 600 700 800 900 t, с

Рис. 2. Графики зависимостей от времени уровня воды в створе плотины (а), объёма (б), расхода (в) и скорости перелива (перехлёста) воды через гребень плотины (г); параметры: H = 100 м, a = 100 м, to = 20 с, vo = 10 м/с, L = 5000 м, xo = 2500 м, Е0 = 500 м

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

а

в

г

Л = 4а = 400 м, Н = 100 м.

Скорость волны, вычисленная по формуле (52), равна с = 25 м/с. Следовательно, численное значение скорости образованной волны точно совпадает со значением, получаемым по общеизвестной классической формуле.

На рис. 2, б представлена зависимость объёма воды, перелитой через гребень плотины, от времени. График имеет ступенчатую форму в связи с тем, что перелив наблюдается в промежутках времени, когда уровень воды у створа плотины превышает уровень у гребня плотины, т. е. при п(0,г) > 0. При п(0,г) ^ 0 волна от плотины полностью отражается и излив воды не происходит.

-4

о 100 200 300 400 500 600 700 800 900 t, с о 100 200 300 400 500 600 700 800 900 t, с

О

О 100 200 0 400 500 600 700 800 900 t, с 0 100 200 300 400 500 600 700 800 900 t, с

Рис. 3. Графики зависимостей от времени уровня воды в створе плотины (а), объёма (б), расхода (в) и скорости перелива (перехлёста) воды через гребень плотины (г); параметры: H =100 м, a = 100 м, to = 20 с, vo = 5 м/с, L = 5000 м, xo = 2500 м, Е0 = 500 м

На рис. 2, в, г приведены зависимости расхода и скорости перелива воды через гребень плотины от времени. Как видно из графиков, с течением времени расход воды уменьшается более интенсивно, чем скорость перелива. Колебание скорости между двумя максимальными (пиковыми) значениями обусловлено колебанием уровня воды у плотины в промежутках времени между двумя максимальными (пиковыми) значениями уровня воды.

Вычислительные эксперименты показали, что количественные значения характеристик образованных волн наиболее чувствительны к скорости вторжения массива в водохранилище. На рис. 3 представлены графики при скорости вторжения в два раза меньшей, чем в предыдущем случае, т.е. v0 = 5 м/с. Приведённые данные указывают на прямо пропорциональную зависимость между амплитудой образованных волн и скоростью вторжения массива в водохранилище.

a б

Г), м w, м3

Результаты настоящей работы получены аналитическими методами решения начально-краевых задач. При этом учитывались непризматическая конфигурация при-плотинной части водохранилища и процесс перелива воды через гребень плотины. Полученные теоретические зависимости могут быть применены для вычисления характеристик образованных волн в узком глубоком водохранилище. Кроме того результаты работы могут быть использованы в качестве тестовой задачи для установления достоверности данных, полученных конечно-разностными методами решения начально-краевых задач.

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

[1] Noda E.K. Water waver generated by lands-lider // J. of the Waterways Harbors and Coastal Eng. Division. 1970. Vol. 96, No. 4. P. 835-855.

[2] Мамрадзе Г.П., Музаев И.Д. Возникновение волн в водохранилище вследствие оползневых явлений // Сообщения АН ГССР. 1971. Т. 64, № 2. С. 115-120.

[3] Гвелесиани Т.Л., ДжинджихАшвили Г.Я. Прогноз максимальной амплитуды плоской волны в водохранилище, возникающей при сейсмогравитационных деформациях // Строительство и архитектура. Сейсмическое строительство. 1980. Вып. 9. С. 27-29.

[4] Остапенко В.В. Численное моделирование волновых течений в Сарезском озере, вызванных катастрофическим обрушением берегового оползня // Вычисл. технологии. 1994. Т. 3, № 8. С. 116-125.

[5] Гогелиани Л., Амброладзе Т. Гидравлическая модель лавинно-оползневого потока, вторгшегося в водохранилище под воздействием сейсмических сил // Гидравлика. Гидротехника. Гидроэкология. Тбилиси: Самигобло, 1996. С. 60-66 (на грузинском языке).

[6] Березин Е.Н., Бейзель С.А. Параллельная реализация алгоритма для расчёта генерации длинных поверхностных волн цунами движением оползня // Вычисл. технологии. 2009. Т. 14, № 1. С. 7-20.

[7] Шокин Ю.И., Чубаров Л.Б. О подходах к численному моделированию оползневого механизма генерации волн цунами // Там же. 2006. Т. 11. Спец. выпуск, посвященный 85-летию со дня рождения академика Н.Н. Яненко. Ч. 2. С. 100-111.

[8] Бейзель С.А., ХАкимзянов Г.С., Чубаров Л.Б. Моделирование поверхностных волн, порождаемых подвижным оползнем, движущимся по пространственно неоднородному склону // Там же. 2010. Т. 15, № 3. С. 35-51.

[9] Шокин Ю.И., Федотова З.И., ХАкимзянов Г.С. и др. Моделирование генерации волн цунами движением оползня с учётом вертикальной структуры течения // Современные методы математического моделирования природных и антропогенных катастроф. Труды VIII Всероссийской конф. Кемерово: Институт угла и углехимии СО РАН, 2005. С. 3-27.

[10] Музаев И.Д., Музаев Н.И. Постановка и решение начально-краевой задачи поверхностных гравитационных волн в водохранилище узкоканьонного типа // Изв. высших учебных заведений. Северо-Кавказский регион. Естественные науки. 2009. № 2. С. 22-24.

[11] Музаев И.Д., Музаев Н.И., Дзебоев Б.А. Математическое моделирование сейсмостойкости плотины с учётом влияния водной среды // Вычисл. технологии. 2012. Т. 17, № 1. С. 90-99.

[12] Сретенский Л.Н. Теория волновых движений жидкости. М.: Наука, 1977. 815 с.

Поступила в 'редакцию 23 октября 2013 г., с доработки — 26 февраля 2014 г.