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

УДК 519.6 ББК 22.11

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

ПРОГРАММНАЯ РЕАЛИЗАЦИЯ

ОДНОМЕРНОЙ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ТРАНСПОРТА НАНОСОВ В ПРИБРЕЖНОЙ ЗОНЕ ВОДОЕМА

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

Ключевые слова: Математическое моделирование, транспорт наносов, аккумуляция и абразия; метод баланса; дискретная модель.

E. A. Protsenko

PROGRAM REALIZATION OF ONE-DIMENSIONAL MODEL OF SEDIMENT TRANSPORT IN COASTAL AREA

Abstract. Problem of building of one-dimensional model of sediment transport in coastal area was considered in this work. Discretization of a mathematical model based on the balance method is made in this work. For the given mathematical model the discrete model is constructed, its program realization is considered, results of numerical experiments are presented.

Key words: mathematical modeling, sediment transportation, accumulation and erosion; a balance method; discrete model.

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

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

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

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

При моделировании процесса перемещения наносов использованы методы вычислительной математики [5, 6], теория разностных схем [4]. При разработке гидродинамической модели движения водной среды использована гидростатическая модель, включающая уравнения движения Навье-Стокса, неразрывности для несжимаемой жидкости, гидростатики [1, 12].

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

Математическая модель транспорта наносов в прибрежной зоне водоема. Исходным уравнением математической модели транспорта наносов является [3, 7, 11]:

= (1)

dt дх

Qs = Amdy/g , (2)

sinS

sin^p " (3)

Ч±--ТЬс

Vs=- ( w '

ть ~ -0,6С11Ь \иь\,

(4)

где е - пористость грунта; Н - глубина дна, отсчитываемая от невозмущенного уровня водоема; <2 - расход наносов; Акр- безразмерные постоянные (в настоящей работе .4=19,5; /?=3); г,, -касательное напряжение на дне; ят - частота волны; ^ - параметр Шильдса (касательное напряжение на дне, записанное в безразмерном виде); с/. р ро - характеристики осадков и воды; С -коэффициент сопротивления твердых частиц; ~иъ - придонная скорость в области неразрушенных волн; х - тангенциальное напряжение на дне водоема; Б - угол, составленный касательной к контуру дна в момент времени /: (р{] - угол откоса грунта в воде; g - ускорение силы тяжести.

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

— х/х—х "■ гДе и, - скорость на критической глу-

/ п 0 ^

половины периода волны: \] =\] х

Ъ Ьс п

бине Н; х - координата верхней границы наката; х - координата точки опрокидывания волн; п

- постоянная (в данной работе п=0,33).

Уравнения (1)-(4) дополняются начальным условием:

Я(х) = Б0х при г = г0.

(5)

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

Н(хн,г) = Нн, хн=#н/я0. (6)

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

Н(хгл > 0 = Нгл ■> Нгл = Л /2-Уравнения (1) - (4) могут быть записаны в компактной дивергентной форме [3, 11]:

(8)

дН д ( дН

---а-

дt дх,

гДе а-к-

Чс

-, Ъ =кть; к =

Атпй

(

1-е

Ти

(А-Ро)я^ •

вт^,, дх

Поставленная начально-краевая задача решается приближенно конечно-разностным методом [4; 6].

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

равномерной сеткой.

со: х{ = х0+1-Ь, \ = 1,2,...,ТУ; ?я=?0 + я-г, п = 0,1,...,М , где и,г - индексы по временной и пространственной координатным направлениям; г, /г - шаги по временной и пространственной координатным направлениям; М, N - количество узлов.

Дискретную модель транспорта наносов запишем в виде:

}{ — ]-["

1+112

■-а

Н"+1-Н,

г-1

Ь" -Ь"

¡+1/2 г-1 / 2

-1/2

(9)

где #г"=#(хгЛ);

, л" — ъп

п" 1+1 /2 г+1 / 2 "

1Ъс

. Т п _ 1 п

■ ' "г+1/2 — К-1+\12ГЬ,1+\12 ■>

к" =

г+1/2

Ашс1

Ти-

т, Н"-Н

ОС__1+1_I

\2

вш ф0

к

1-е

Ьс

3

Таким образом, построена двухслойная неявная конечно-разностная схема. Математическая модель гидродинамики. В рамках теории мелкой воды исходными уравнениями движения водной среды являются [1; 8; 11; 12]:

- система уравнений Навье-Стокса:

1 ' '

и[+ии'+м>и'=--Р' + ии' + пи' , (10)

I X 2 X ' X V > 2 ?

р

1 '

+ ин>'х + УУУУ'2 ---Р' + /ш' + Щ>\ >

р х 1

- уравнение неразрывности для несжимаемой жидкости:

и'+м>'= О, (11)

- уравнение гидростатики:

Р = Рё г + ^

(12)

где ^(х,/) - отклонение свободной поверхности жидкости от равновесного состояния; Р - гидростатическое давление; р - плотность жидкости; //. // - составляющие коэффициента турбулентного обмена; и, ^ - компоненты вектора скорости; х - горизонтальная координата; / - время; g -

ускорение свободного падения.

Система уравнений (10) - (12) дополнена следующими граничными условиями:

- на дне задается условие непроницаемости и трения: руГ!и'п — тъ (/), м? = 0,

- на поверхности задается подъем уровня и ветровые напряжения:

- на боковых границах задается условие свободного выхода: и'п = 0, СИ = 0 •

При помощи интегро-интерполяционного метода получена аппроксимация математической модели гидродинамики (10) - (12), которую можно представить следующими уравнениями:

- для расчета поля скорости на промежуточном временном слое:

Л. - и.

2 1 Ьг+1/2 7+1/2

И.

и1+1-иг

2/г

+

Т-Г £ ^г и £

+ -"г-1/2 + Ьг-1/2 Мг-1/2 = -"й-1/2 + £г

2/2

г+1/2 1 Ьг+1/2 М+1/2

Щ+х-Щ

-И П I Г

Л г-1/2 Ьг-1/2 М-1/2 ^2

- для расчета функции возвышения уровня:

р р

(

н1

^г+1/2 + ¿¿+1/2

^г+1/2 + ¿¿+1/2 "¿+1/2

^ _ тт . ^ ^-1

77 г-1/2 Ьг-1/2 ,_2

/г2

М,-,, „ Н, | 2 + Мг-1/2 .

А

(13)

(14)

- для уточнения поля скорости на следующем временном слое:

Н,+ ^ ^-

тт , £ ^+1 ^ . Т_Г . £ ь!

¿+1/2 ^¿+1/2 „, ЛН/21"ЬН/2

(15)

2/г "" 2/г

\ /

Здесь и - значение поля скорости на текущем временном слое; и - значение поля скорости на предыдущем временном слое; й - значение поля скорости на промежуточном временном слое й -ай+ \ - с и, <У - вес схемы [4, 9].

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

метод вариационного типа [6; 7; 11; 12], который рекомендуется применять при решении задач транспорта наносов большей размерности.

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

Областью применения данного программного комплекса является построение турбулентных потоков несжимаемого поля скоростей водной среды и прогнозирование возможных сценариев изменения геометрии донной области мелководных водоемов.

Для функционирования программы на вычислительной платформе необходима предварительно установленная операционная система не ниже Windows XP. Программа написана на языке C++.

Разработанный комплекс программ состоит из четырех программных компонент, предназначенных для:

- расчета поля скорости на основе волнового уравнения;

- расчета поля скорости на основе волнового уравнения с учетом нелинейных процессов;

- расчета поля скорости на уравнения Навье-Стокса;

- расчета изменения геометрии дна водоема.

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

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

- блок построения сеточных уравнений для волнового уравнения в соответствии с трехслойной конечно-разностной схемой;

- блок построения сеточных уравнений для волнового уравнения с учетом нелинейных процессов;

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

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

- блок расчета поля скорости с учетом возвышения уровня (результатом работы данного блока является расчет значений поля вектора скоростей на следующем временном слое);

- блок расчета сеточных уравнений прогонкой;

- блок вывода значений поля скоростей в файлы.

В программном комплексе задаются следующие параметры:

- размерность сетки по пространственной координате;

- шаг по времени и по пространственной координате;

- расчетный отрезок времени;

- частоту и амплитуду волновых колебаний;

- критическое значение касательного напряжения, при котором начинается перемещение наносов.

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

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

На рисунках 1 - 5 представлена динамика изменения рельефа дна и функции возвышения уровня в случае наклонного дна. Максимальная глубина водоема составляла 1 метр, минимальная - 50 см.

0.21-

0.14 0.08 0.02

- 0.04 - 0.1

- 0.16

- 0.22

- 0.28

- 0.J4

j - 0.4

- 0.46

- 0.52

- 0.58

- 0.61 - 0.7

- 0.76

- 0.82

- 0.88-

- 0.91-

- I1-

0 400 800 1.2xl031.6vl032*103 2.4vl0%.Sxl0^.2vl0^.6>,1034xl03 4.4>,1034.Sxl0;i.2xl035.6xl03 6<103 6.4х1036.8>1037.2х1037.6>103 8>103

Рис 1. Функция возвышения уровня через время равное половине периода волны в случае динамически изменяемого рельефа дна

0.21-

0.14 0.08 0.02

- 0.04

- 0.1 —

- 0.16 —

- 0.22-

- 0.28-

- 0.34-

) - 0.4

- 0.46-

- 0.52-

- 0.58-

- 0.64-

- 0.7-

- 0.76-

- 0.82-

- 0.88-

- 0.94-

- I1-

0 400 800 1.2а1031.6х1032к103 2.4*10^.8*10^.2*1033.6*1034>.1034.4х1034.8x10%.2а1035.6а103 6х103 6.4х1036.8*1037.2*1037.6*103 8х103

Рис. 2. Функция возвышения уровня через время равное периоду волны

0.21-

0.14 0.08 0.02

- 0.04 -0.1 —

- 0.16 — -0.22 — - 0.28 —

- 0.34 — ¡-0.4

- 0.46 —

- 0.52 —

- 0.58 —

- 0.64 — -0.7-

- 0.76 — -0.82 —

- 0.88 — - 0.94 —

- I1-

О 400 800 1.2а1031.6х1032х103 2.4*10^.8*10%.2*1033.6*1034х1034.4х1034.8х10%.2А1035.6А103 6х103 6.4х10;^.8*1037.2*1037.6*103 8А103

Рис. 3. Функция возвышения уровня через полтора периода волны

0.21-

0.14 0.08 0.02

- 0.04 -0.1 —

- 0.16 — -0.22 — - 0.28 —

- 0.34 — ¡-0.4

- 0.46 —

- 0.52 —

- 0.58 —

- 0.64 — -0.7

- 0.76 — -0.82 —

- 0.88 —

- 0.94 —

- I1-

О 400 800 1.2х1031.6х1032х103 2.4*10^.8*10%.2*1033.6*1034х1034.4х1034.8х10%.2х1035.6х103 6х103 6.4х10^.8*1037.2*1037.6*103 8х103

Рис. 4. Функция возвышения уровня через два периода волны

0.21-

0.14 0.08 0.02 -0.04 -0.1 — -0.16 — -0.22 -0.28 — -0.34 — ¡-0.4 -0.46 —

- 0.52 — -0.58 — -0.64 —

-0.7

- 0.76 — -0.82 — -0.88 — -0.94 —

- I1-

О 400 800 1.2*10\бх1032х103 2.4х1032.8*10^.2х1033.6х1034х103 4.4*1034.8х10%.2х1035.6х103 6х103 6.4х1036.8х1037.2<1037.6*103 8<103

Рис 5. Функция возвышения уровня через два с половиной периода волны

На рисунках 6 - 10 представлена динамика изменения рельефа дна в случае наклонного дна.

0.014 _ , „-3

2.4x10 3.2x10

4.8x10 5.6x10

Рис. 6. Изменение рельефа дна через время равное половине периода волны

0.051-

0.041 0.032 0.023 0.014 5x10"3

-4x10"3----------

- 0.013 -0.022 -0.031

_от1----------

О 800 1.6x103 2.4x1 о3 3.2x103 4х103 4.8х103 5.6х103 6.4x103 7.2х103 8х103

Рис. 7. Изменение рельефа дна через время равное периоду волны

0.051-

0.04 0.03 0.02 0.01 ] 0

- 0.01 - 0.02

- 0.03

- 0.04

- 0.05'----------

О 800 1.6х103 2.4х103 З.2х103 4х103 4.8х103 5.6х103 6.4х103 7.2х103 8х103

Рис. 8. Изменение рельефа дна через полтора периода волны

0.051-

0.04 0.03 0.02 0.01 0 0.01 0.02 0.03 0.04

0.05----------

О 800 1.6х103 2.4x103 З.2х103 4х103 4.8х103 5.6х103 6.4х103 7.2х103 8х103

Рис. 9. Изменение рельефа дна через два периода волны

0.05 0.04 0.03 0.02 0.01 0 0.01 0.02 0.03

0 800 6.1" 1I I" 3.2/10"' I-1" 4.8'103 5.6x10"' 6 1-1" 7.2Л03 И. I"

Рис. 10. Изменение рельефа дна через два с половиной периода волны

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

- транспорт наносов осуществляется в направлении распространения волн, как в случае плоского дна, так и в случае наклонного;

- приближенная математическая модель формирования наносов в прибрежной зоне водоема с определенной точностью совпадает с аналитическим решением [2].

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

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

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

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Сравнительный анализ классических и неклассических моделей гидродинамики водоемов с турбулентным обменом / Е. В. Алексеенко и др. // Известия ЮФУ. Технические науки. -№ 8(97). - Таганрог: Изд-во ТТИ ЮФУ, 2009. - С. 6-18.

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

3. Проценко, Е. А. Модель и алгоритмы решения задачи о транспорте наносов // Известия ЮФУ. Технические науки. - № 8(97). - Таганрог: Изд-во ТТИ ЮФУ, 2009. - С. 71-75.

4. Самарский, А. А. Теория разностных схем: учебное пособие для вузов / А. А. Самарский. - 3-е изд., испр. - М.: Наука. Физматлит, 1989. - 616 с.

5. Самарский, А. А. Численные методы решения задач конвекции-диффузии / А. А. Самарский, П. Н. Вабищевич. - 4-е изд. - М.: Либроком, 2009. - 248 с.

6. Самарский, А. А. Численные методы математической физики : учеб. пособие / А. В. Гулин, А. А. Самарский. - 2-е изд. - М.: Научный мир, 2003. - 316 с.

7. Сухинов, А. И. Построение и исследование дискретной математической модели расчета прибрежных волновых процессов / А. И. Сухинов, Е. Ф. Тимофеева, А. Е. Чистяков // Известия ЮФУ. Технические науки. - № 8(121). - Таганрог: Изд-во ТТИ ЮФУ, 2011. - С. 22-32.

8. Сухинов, А. И. Численная реализация трехмерной модели гидродинамики для мелководных водоемов на супервычислительной системе / А. И. Сухинов, А. Е. Чистяков, Е. В. Алексеенко // Математическое моделирование. - 2011. - Т. 23. - № 3. - С. 3-21.

9. Сухинов, А. И. Оценка погрешности решения уравнения диффузии на основе схем с весами /

А. И. Сухинов, А. Е. Чистяков, Ю. С. Бондаренко // Известия ЮФУ. Технические науки. -№ 8(121). - Таганрог: Изд-во ТТИ ЮФУ, 2011. - С. 6-13.

10. Сухинов, А. И. Построение дискретной двумерной математической модели транспорта наносов / А. И. Сухинов, А. Е. Чистяков, Е. А. Проценко // Известия ЮФУ. Технические науки. -№ 8(121). - Таганрог: Изд-во ТТИ ЮФУ, 2011. - С. 32-44.

11. Сухинов, А. И. А. Двумерная гидродинамическая модель, учитывающая динамическое перестроение геометрии дна мелководных водоемов / А. И. Сухинов, А. Е. Чистяков, Е. А. Проценко // Известия ЮФУ. Технические науки. - № 8(121). - Таганрог: Изд-во ТТИ ЮФУ, 2011. - С. 159-167.

12. Чистяков, А. Е. Трехмерная модель движения водной среды в Азовском море с учетом транспорта солей и тепла // Известия ЮФУ. Технические науки - № 8(97). - Таганрог: Изд-во ТТИ ЮФУ, 2009. - С. 75-82.

УДК 531.38 ББК 22.21

С. А. Шретер

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ПЛОСКОПАРАЛЛЕЛЬНОГО ОБТЕКАНИЯ ПЛАСТИНКИ ПОТОКОМ ВОЗДУХА1

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

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

S. A. Shreter

MATHEMATICAL MODEL OF PLANE-PARALLEL FLOW PAST A PLATE AIR FLOW

Abstract. Consider the problem of the behavior of hybrid system in a stream of air. The idea of the solution consists in reducing the original equilibrium equation of Kirchhoff to the system of equations Hamiltonian type, where the Hamiltonian function is normalized in a certain number of members.

Key words: Hamiltonian approach, Birkhoff transformation, bending the rod, the mathematical model, the aerodynamic forces.

Постановка задачи. Рассмотрим эксперимент по определению зависимости аэродинамических параметров от конечных перемещений точек упругого стержня. Стержень помещают в набегающий поток воздуха. Нижний конец стержня жестко защемлен, к его верхнему концу жестко прикреплена абсолютно твердая пластинка. Предполагается, что поток воздействует только на пластинку, изгиб стержня происходит в одной плоскости. Начальное положение стержня определяется заданием угла наклона стержня в=ц/ по отношению к скорости набегающего потока (см. рис.1). Пластинка (также как и тело произвольной формы) отличается от точки тем, что направление силы воздействия потока на тело R может не совпадать с направлением вектора V скорости потока. Поэтому силу воздействия R потока на пластинку представим в виде суммы двух векторов: R=S+P, где S - сила сопротивления, Р -подъемная сила, SUF, PJJ7. Для аэродинамических сил возьмем зависимости [7]:

S=s{a)VV% . /' /.u/iCt/dV', ,

где р - плотность воздуха, а - угол атаки пластинки, i ( Ш') - единичный вектор, лежащий в плоскости пластинки. Функции s(a), р(а)~ коэффициенты аэродинамических сил, зависят от формы и размеров пластинки и определяются экспериментально. Общие свойства этих функций:

s(a) > sQ >0,а р(а) меняет свой знак при а=0 и а=ж/ 2 (р(Р)=р(х/2)=0).

1 Данная работа выполнена при финансовой поддержке государственного задания Министерства образования и науки РФ

ФГБОУ ВПО «ТГПИ имени А.П. Чехова» по проекту № 1.1885.2011, тема: «Математическое моделирование статики и

динамики гибридных механических систем и идентификация их параметров», научный руководитель А. А. Илюхин.