Двумерная конечно-разностная модель формирования наносов в прибрежной зоне водоема и ее программная реализация Текст научной статьи по специальности «Математика»

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

Е.А. Проценко

ГОУ ВПО ТГПИ, Таганрог

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

Проблемам математического моделирования береговых процессов посвящены работы Э.Б.Бертмана, И.В.Попова, Н.В.Пыхова, И. О. Леонтьева, К.В.Гришанина, И.Г.Кантаржи, А.В.Караушева, Ч.Т.Янга и мн.др.

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

В данной работе освещены вопросы построения пространственно-двумерной модели транспорта наносов. Описанная математическая модель используется для численного моделирования динамики аккумулятивного берега. В соответствии с критерием «крутизны» используется уравнение непрерывной модели формирования наносов.

Уравнения процесса перемещения наносов записываются в виде:

где Н - глубина дна, отсчитываемая от невозмущенного уровня водоема; 8 -пористость грунта; Q - расход наносов; х, у - направления дифференцирования; тъ -

касательное напряжение на дне; ТЪс - критическое значение касательного напряжения, при котором начинается перемещение наносов, А ив - безразмерные постоянные (в работе А = 19,5; в=3), - частота волны, ё- характеристика осадков.

Формулы для расхода наносов:

где /0 - параметр Шильдса (касательное напряжение на дне, в безразмерном виде); р1, р0 - плотности твердых частиц и воды; С - коэффициент сопротивления твердых

частиц; иъ - придонная скорость в области неразрушенных волн. Значения иъ и тъ

интерпретируются как усредненные в течение половины периода волны. Трансформация параметров волн до глубины опрокидывания рассчитываются по линейной теории волн, в области прибойного потока - по эмпирическим зависимостям.

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

д дх ду

Ub = Uc

b bc

x - xс

V n су

(5)

где Ubc - скорость на критической глубине Н; Xn - координата верхней границы наката; х0 - координата точки опрокидывания волн; п - постоянная (п=0,33).

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

берега. Предполагается, что в начальный момент времени (ї=0) к откосу, сложенному из осадочных пород, подходят монохроматические волны. Ось ОН, представляющая изменение глубины наносов, отсчитывается от невозмущенной поверхности жидкости и направлена в сторону углубления. Схематически данная система координат представлена на рисунке 1:

Изначально считаются заданными следующие параметры: Б - начальный уклон дна; И0, А0 - высота и длина волны; т - период волны; й, р15 р0 - характеристики осадков и воды; Т - длительность шторма.

Система уравнений для параметра Шильдса принимает вид: = /0

Г \

sin S 1 ±-------

V &Рс J

(6)

где /Б - параметр Шильдса для наклонного дна; Б(х,у,/) - угол, составленный касательной к контуру дна в момент времени /; (ро - угол естественного откоса грунта в воде.

При прочих равных условиях транспортирующая способность потока, переносящего наносы вверх по откосу, меньше, чем транспортирующая способность потока над горизонтальным дном и параметр Шильдса записывается со знаком «-». Соответственно

расход потока, переносящего частицы в направлении к берегу: q

= Am d/

1 -

sin S

tgp

(7)

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

Подставим уравнение (7) в исходное уравнение наносов (1):

(1 -s)

dH

~dt

д_

dx

Arnd

1-

sinS

tgp_

д_

ду

Arnd

1-

sinS

tgp

с

= с.

При решении пространственной задачи о переформировании берегов полагаем: sin S - S

arctg.

дH

дx

+

ґд^2

ду

дH

дx

+

ґд^2

ду

(8)

Графически данный вывод представлен на рисунке 2.

дН

ну

Подставим (9) в уравнение (8), имеем:

л ^дH д

(1 -є)------+ —

v 7 ді дх

(

ЛтсС

1-

8ІИ Б

в

У

( (

д_

ду

(

ЛтсС

1-

8ІИ Б

в

л ЄдН д

(1 -є)------+ —

ді дх

ЛтсС

1-

дН У (сН_ дх ) I ду

^Фо

ду

ЛтсС

1-

дН У2 (дНЛ 2

дх ) I ду

в

^Фо

= 0.

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

л Ч дН д

(1 -є)-------+ —

ді дх

ЛтсС

дН У (дН42 +

1 уу дх ) У ду ^Фо

ду

ЛтсС

дН У (дН42

1 уУ дх ) у ду

^Фо

= о,

(1 -є) дН + (А + А

ді У дх ду

ЛтсС

.9 / \2

дН У (дН

1-

дх ) У ду

&Фо

V

(

= о,

л ч дН ( д д

(1 -є)-------------+ 1 — + —

ді У дх ду

ЛтсС

1-

дН У2 +(Н2 дх ) І ду

Ув-1(

іЯФо

V

л ^дH ( д д

(1 -є)------+ 1 — + —

ді У дх ду

1-

ЛтсС

1-

дН У2 + ( дН дх ) У ду

tgФо

дН У2 | ( дН дх ) У ду

У V -\Р-1 г

1-

іЯФо

дН У2 + ( дН дх ) У ду

= о,

дН У2 (дН42

аТ і +Уа^1 'Ф

= о.

л л дН д

(1 -є)-------+ —

ді дх

Л тС

1-

дН У2 (дН42 +

дх ) У ду

V-

&Фо

1-

дН

дх

дН У2 (дН42 , +

дх

ду

і8Фо

д

+

2

2

+

дх

ду

'ди Л2 (ди Л2

+1 I

кдх ) I ду )

1 -■

дИ

ду

дИ Л ( дИ ,

,аГ МзГ 1

ду

Продолжая преобразования, получим:

= 0.

п ) дИ 8

(1 - є)--------------------

д, дх

Атеі

і Л

дИ_

дх

2 ( дИ Л

Лв-1(

ду

ї?,Фо

дИ

дх

дИ

дх

2 І дИ Л

+ -------

ду

&00

_д_

ду

Аті

1 л

дИ

дх

----------------\в-і /

2 / \2 '

2 ! яи \

дИ

ду

1 --

( дИ Л ду

дИ

дх

2 ( дИ Л

( ( (дИЛ2 +|'8Я 12 ' в-1 Л ( (

д д х Ат сі 1 -^ 1 дх ) У ду ) д ду Аті

о

V V ) ) V 1у~

ду

і-V

,ё^о

дИ Л2 ( дИ 42

\ +

д х

ду

\Р-1\

= 0.

(10)

Представим уравнение (10) в сокращенной, удобной для восприятия форме:

аь є >

2 Л д( (И л2Л

где

Сх дх

+ с, + су = 0

Ь =-

(11)

1

' дИ

2 ( дИ Л

с =-

дх ( ( Аті

ду

(12)

ду

(1 -є)

1 л

дИ 'Л2 І дИ Л дх \ + V ду

■ Лв-1 Л

^0

д

д

+

2

2

+

2

2

Итак, уравнение (11) описывает переформирование дна в случае, если все частицы наносов двигаются в сторону берега.

Уравнение неразрывности дополняется начальным условием:

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

Таким образом, имеем непрерывную двумерную математическую модель формирования наносов в прибрежной зоне водоема (11) - (15).

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

Предполагаем, что используется временная сетка с постоянным временным шагом

Границы области О определяются, исходя из физических соображений. Одна из границ, близкая к береговой линии, но не совпадающая с ней, определяется как верхняя граница наката, где средняя скорость волн обращается в ноль, и берег не подвергается деформациям: Н(хн,у,0 = Ин, хн = Ин /£0.

Другая (жидкая) граница определяется как условная граница «глубокой воды», где влияние волнения на подъем взвешенного вещества пренебрежимо мало:

Другие граничные поверхности области G определяются пользователем. Для них предполагается, что известны потоки наносов в направлении к нормали к ним.

Будем вначале для простоты считать, что область G прямоугольник - покрыта равномерной сеткой Шк: (ок =ах хау, где

(13)

Н (хн, У,1) Нн, хн Нн / V

На границе «глубокой воды» глубина также не изменяется:

н (хгл,у,1) = н гл , н гл =Л/2.

(14)

(15)

а>Т = {іп = пт,0 < п < К}.

н(хГЛ,у, 0 = нГЛ , нГЛ =Л)/2.

Т х

(16)

где а =

V

У

V

У

1

1

2

а1-1х (х, У}) = а1 (х - К, у}), а2-1у = а2 (хі, У} - Ку ) ,

ъ1-1х (хі, У; ) = Ъ1 (хі - К, у}-), ъ2-1у (хі, У}- ) = Ъ2 (хі, У}- - КУ )

с = — х 2

Люё

1

с = — у 2

(1 -е)

( ( Лшё

1Л

*ЯЯ>о

Люё

(1 -е)

1 л

(Н )2 + ( Н

у

(1 -е)

1-

(НУ )2+(Н

-у-1 л

1

+—

2

( ( Л<а(1

(1-е)

1 -

(НУ) + (Н

(17)

;у

В соотношениях (16) и (17) используются следующие обозначения [2]: /х = I, (, У-)-^; /- = /- (, у.)

/ (х, + К У')-/(,, У,) ,_,(„ .Л- 1 (, • У>)-1 (, - К' У>-

К

/о = /о (, , У, -

2к

; /у = /у (, у, -

/ ( ,, у, + ку )-1 ( ,, у, ),

= ,(х у)-1 (х,у,)-1 (х',у,-ку). г = г (х у )-1 (,■,у, + К)-1 (х,у,-К)

1у 1у\ ',У]) к ’ у °у ', ]' 2к ’

символ " /" означает, что значение сеточной функции / рассматривается на текущем -верхнем временном слое (номера п +1); функции без символов над ними предполагаются заданными на предыдущем временном слое (номера п ).

Такая дискретизация позволяет получить разностную схему, имеющую погрешность аппроксимации О (2 + ку 2 + т).

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

Расчет изменения глубины на следующем временном слое производится по следующему алгоритму: вычисление коэффициентов сеточного уравнения; определение начальной невязки и ее нормы; итерация методом минимальных невязок или попеременнотреугольным методом; вычисление невязки; условие выхода из цикла по заданному параметру е.

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

у.к+1 ,к

Канонический вид уравнения: х - х + л,к = /к = 01 (18)

Значение невязки: гк = / - Лхк, к - шаг по времени.

Умножим обе части (18) на матрицу А справа: —Л,~+ЛЛХк = Л/'

Тк+1 ‘

(19)

Меняя знаки и группируя слагаемые соответствующим образом, продолжаем преобразования: (-(х + /)~{-Лх +1) + Л[-Л,к + /) + Л/ = Л/ и, учитывая равенство

х

у

(19), имеем:

rk+\ _ rk

+ Ar = 0.

(20)

k+1

Итерационный параметр Tk+1

,k+1 к a к

к+1

, будем выбирать из условия минимума невязки r по

норме r

r

rk+1

тк+1 ■Ar

^ тМ =((rk _Ти+М ),(rk _4+1AjJk)) =

/|2 _ 2zk+1(Ark, rk) + (tJ2| Ark Продифференцируем <р(тк+1) по тк+1, имеем:

_2(Ark, rk) + 2Tk,1||Ark|f = 0, . =(Ar‘, r‘)_(Ark, rk)

- ||Агк| |2 (гк, Агк)

Найдем вторую производную: <р"(тк+1) = 2||Агк||2 > 0.

Следовательно, Тк+1, удовлетворяющее равенству, - это точка минимума. Оценим скорость сходимости метода минимальных невязок.

к+1

_Т+1Ar

<

_Т Аг

<|\Е _То A|| ’I г

’ Т0

Y1 +Y2

-, гЕ < A < Y2Е.

Предположим, что у2>>у1, тогда т0 «—; y2

Y2

IIе _т 4

Таким образом,

Е _—A Е _ — Е 1

Y2 IAII 2 = 2

к+1

2 <fi

В ходе выполнения данной работы был разработан программный комплекс, реализующий ввод исходных данных и численное решение разностных уравнений методом минимальных. Ввод в программном комплексе реализован в виде трех блоков. Первый блок - блок заполнения параметров предназначен для ввода параметров волны, а именно, длины и амплитуды волны. Следующий блок предназначен для заполнения первоначальных параметров дна (первоначального угла наклона) и пористости грунта. Завершающий блок позволяет определить геометрию области, в том числе параметры сетки. Также необходимо указать длительность действия ветра с соответствующими значениями параметров волны, то есть интересующий промежуток времени. Если значения полей длины и высоты волны не соответствуют «критерию крутизны», то есть процесс относится к деструктивному режиму (абразия), то при нажатии на кнопку «Start» появится соответствующее сообщение. Аналогичное сообщение появится, если одно из полей не заполнено.

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

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

Были проведены исследования для различных длин волн, наклонов дна и пористости грунта и с погрешностью 10-25 % данные численных экспериментов совпали с существующими данными образования наносов в восточной части Таганрогского залива.

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

Propates о! a wave Wave-length Surge height 9 і

Paosily 0.1

Angle of dip 15

N1 20

N2 20

Stum time [То

Рис. 3. Примеры работы программы для различных значений начальных параметров

О 5 10 15

С

Рис. 5. Распределение глубин для заданной области

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

зо

Результаты данной работы могут быть применены на начальных этапах прогнозирования формирования прибрежной зоны водоемов.

Литература

1. Леонтьев И.О. Прибрежная динамика: волны, течения, потоки наносов. - М.: Геос., 2001, -272 с.

2. Самарский А.А., П.Н. Вабищевич. Численные методы решения задач конвекции-диффузии. - М.: Едиториал УРСС, 1999. -248 с.

3. А.И. Сухинов, Е. А. Проценко. Модифицированный попеременно-треугольный метод решения разностной краевой задачи Дирихле для уравнения эллиптического типа в прямоугольнике с линейной функцией источника. Методы и алгоритмы прикладной математики в технике, медицине и экономике: Материалы VII Междунар. нау. - практ. конф., Новочеркасск, 2 фев. 2007 - Новочеркасск: ЮРГТУ, 2007. - Ч. 1. С.6-8.