Численная модель взвесенесущего потока для Новосибирского водохранилища Текст научной статьи по специальности «Физика»

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

Том 15, № 2, 2010

Численная модель взвесенесущего потока для Новосибирского водохранилища*

В. А. Шлычков

Институт водных и экологических проблем СО РАН, Новосибирск, Россия

e-mail: [email protected]

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

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

Введение

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

Процесс заиления водохранилищ в натурных условиях весьма сложен и в значительной степени индивидуален для каждого водоема. Закономерности перераспределения взвешенных наносов в водоеме зависят от целого ряда факторов: 1 — скорости седиментации частиц (т. е. фракционного состава наносов), 2 — стратификации и интенсивности турбулентного перемешивания, 3 — локальной структуры продольных скоростей, 4 — наличия упорядоченных вертикальных токов, связанных с обтеканием неровностей дна, 5 — морфометрических параметров — длины, конфигурации водного зеркала, размеров поперечников и др.

Надежные и универсальные методы прогноза деформаций ДОННОГО профиля В НЭ/~ стоящее время отсутствуют. Существующие инженерные схемы расчета [?] имеют эмпирический характер, а применяемые в практике балансовые методы основаны на анализе натурных данных и дают лишь интегральные характеристики объемов заиления. Вместе с тем тонкие процессы переноса полидисперсных взвесей могут быть описаны с помощью средств численного моделирования динамики водотоков в сочетании с алгоритмами переноса тяжелой взвеси.

* Работа выполнена в рамках интеграционного проекта СО РАН № 23-09.

© ИВТ СО РАН, 2010.

Применение трехмерной математической модели для описания названных процессов применительно к Новосибирскому водохранилищу, по-видимому, не является оптимальным. Этот водный объект при длине около 180 км (по судовому ходу) имеет сравнительно небольшую среднюю ширину 3.8 км, т. е. в большей своей части носит выраженный продольно-русловый характер без заметного меандрирования и широких пойменных заливов (за исключением Бердского залива в приплотинном сегменте водохранилища). Осредненные по вертикали двумерные (так называемые плановые) модели, широко используемые для расчета динамики влекомых наносов и прогноза русловых деформаций в реках, малопригодны при оценке заиления водохранилищ, так как в последнем случае ведущими являются механизмы переноса взвешенных фракций, и поэтому важна вертикальная детализация процессов массообмена. По этим причинам в качестве рабочего инструмента принята двумерная продольно-вертикальная модель, полученная из уравнений гидродинамики осреднением по поперечному к водотоку направлению. Отметим, что близкий подход использовался для расчета переноса взвеси в эстуарии [?].

Цель Н clC ТОЯ ТД6 и работы состоит в формулировании математической модели течений на основе уравнений гидротермодинамики с поперечным осреднением и получении оценок скорости заиления дна в Новосибирском водохранилище.

1. Модель динамики водохранилища

Пусть у = Ь(х, г) — ширина руслового потока, ось х направлена горизонтально вдоль

основного течения, ось г — вертикально вверх. Осредненные по поперечнику (в на-

у

продольно-вертикальной плоскости будут иметь вид [?]

С

дЬи дЬии дЬит ,дС аЬ д [ . , д , ди д , ди ^

"аТ" + ~я— + ^— = -^^Г + ЬКХ^~ + ЬКХ — - РЬ,

дъ дх дг дх рдх] дх дх дг дг

X

дЬи + дЬи ...

дх дг

где и,т — горизонтальный и вертикальный компоненты вектора скорости; г = £(х, Ъ) — уравнение свободной поверхности; р — средняя плотность воды; р — возмущение плотности; Гь — сила трения вдоль бортов русла; Кх, Кх — коэффициенты горизонтального и вертикального турбулентного обмена; — боковой приток массы.

Обозначим г = гь(х) — уравнение границы дна и рассмотрим краевые условия. На линии дна зададим условия

ът ди II дгь

Кх —— = оа \и\и, и = —— и при г = гь, (2)

д г д х

где оа — коэффициент сопротивления. На свободной поверхности краевые условия имеют вид

ди Тх д( д( . ,

Кх —— = —, ——+ и— = и при г = С, (3)

д г р дЪ дх

Тх

Расчет вертикального турбулентного обмена проводится на основе уравнения баланса кинетической энергии турбулентности е и уравнения для скорости ее диссипации е, которые согласно [?] имеют вид

дЬе дЬие дЬте I д , де д , де . , т "аТ + ~я— + = ЬКХ1Г + ЬКХ — | + ЬКХ3 — Ье,

дъ дх дх \ дх дх дх дх

дЬе дЬие дЬwе ( д 1 де д 1 ^ де \ 1 е т ,е'

"тгг + --+ —^— = о-в -7— ЬКХ— + — ЬКг— + еф-Кг3 - оф—,

дъ дх дх \ дх дх дх д х е е

К7

е е

дТ

(4)

где ае, оо2, о3 — эмпирические постоянные; 3 = щ — двт ^--скорость пополнения

энергии турбулентности, Т — отклонение температуры от среднего профиля, — коэффициент термического расширения воды. Последнее соотношение в (??) выражает гипотезу замыкания Колмогорова.

Для уравнений (??) зададим следующие граничные условия:

д е е3/2

т— = 0, е = ех- при х = хЬ,

дх хзЬ

де

аеКх— = о8 дх

ТХ

рт

е3/2

е = ох- при х = (,

хз(

(5)

где хзЬ, — шероховатости дна и свободной поверхности, ох,о3 — эмпирические посто-

При задании краевых условий на верхнем по потоку сечении х = х1 предполагается отсутствие значимой неоднородности течения (например, за счет резкого изменения ДОННОГО рельефа у входного створа). Вследствие этого в уравнениях (??), (??) можно пренебречь горизонтальными вариациями искомых полей и рассмотреть следующую

х

Ьд С

Ь7Г

дх

д , ^ д и +—ЬКХ—

д х д х

0,

Х = Х1

д д е

—ЬКг — + ЬКХ 3 — Ье

д х дх

д д е е е2

0, —ЬК2— + о2Ь-К, 3 — о3Ь— дх д х е е

0.

(6)

Для уравнений (??) ставятся краевые условия из (??), (??), (??)• Полученные по одномерной модели функции и1,е1,е1 служат граничными значениями в краевом условии на боковом сечении 1:

и = и1, е = е1, е = е1 при х = х1. На выходном (вниз по потоку) сечении х = х2 поставим условия

ди де де

7— = 0, 7— = 0,7— = 0 пр И х = х2.

дх д х д х

(7)

Горизонтальное перемешивание описывается с помощью модели Смагоринского в реализации [?].

Система уравнений для определения геометрии свободной поверхности формиру-

г

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

хода Qi = ЬгкАг, а из уравнения неразрывности следует нестационарное СОСТАВ

ношение, связывающее Qi и Полученная система относительно Q,( является ги-

дЪ

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

Q = Ql при х = х1, ( = при х = х2, (9)

где Q1 — расход воды на границе, (2 — заданная отметка уровня на выходном створе.

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

Рассмотрим способ ЗЭДсШИЯ коэффициента донного трения в (??)■ При решении задач о движении наносов адекватный расчет напряжения трения является важным для определения превышения сдвиговых напряжений над критическими, при которых развиваются явления взмыва. В литературе используются разные гипотезы для расчета этой величины. Так, в экспериментах Никурадзе для условий с полным проявлением шероховатости получено оа = 0.014. В работе [?] принято значение оа = 0.14. Однако расчеты показывают, что коэффициент трения при квадратичном режиме сопротивления не является константой, а зависит от параметров течения — его глубины и скорости, а также от шероховатости дна.

Задача определения оа решалась численно на основе интегрирования уравнений (??) и сопоставления полученных значений придонной скорости с теоретическим значением, следующим из логарифмического профиля скорости [?]

1.25и

иь = -

*| 6а5А

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

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

где ¿ь — средний диаметр донных отложений, а,Ь,с — постоянные, определяемые методом оптимизации. Приведенная формула использовалась при реализации краевого условия (??).

Для проведения расчетов рассмотрим участок Новосибирского водохранилища длиной 130 км, примыкающий к плотине ГЭС. Расчет течения выполнен для расхода Ql = 3000 м3/с, соответствующего периоду весеннего паводка, и при устойчивой стратификации. На рис. ??, а показаны расчетная геометрия свободной поверхности и установившееся распределение продольной скорости вдоль основного водотока. Вниз по течению ширина водного зеркала увеличивается, в силу чего в приплотинной части размер поперечника достигает 20 км. По этой причине скорость течения падает, уменьшаясь от 0.8 м/с вблизи входного створа (х = 0) до 0.1 м/с на озерном участке водохранилища. Рисунок 1, б иллюстрирует структуру вертикального коэффициента турбулентного обмена. Зона максимальной интенсивности турбулентного перемешивания располагается в средних слоях воды и уменьшается с приближением к донной и свободной поверхностям, что соответствует известным теоретическим положениям и расчетам.

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

а

20

10

15

О

5

Ш

40

60

8©

100:

.V, км

Рис. 1. Расчетные поля продольной скорости и, м/с (а) и коэффициента вертикальной турбулентности Кг, см2/с (б) в Новосибирском водохранилище. Горизонтальными стрелками показаны профили скорости па отдельных вертикалях

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

К 6

шероховатости Маннинга п так С = —. Коэффициент Шези отражает суммарное со-

п

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

Вертикальная детализация течения позволяет в явном виде рассчитать суммарный

дЬ\и\и

по х вклад сил трения и путем сопоставления этих значений с агрегатом ———, фигу-

С2

рирующим в уравнениях Сен-Венана, оценить параметры С8и п. Таким образом, для определения С используется уравнение

/ д 1 ^ д и д 1 ^ д и (—ЬКх+ 1- ЬКХ— 1 _ \дх дх дх дх

С2 дЬ\и\и '

где и — средняя по х продольная скорость, вычисляемая по расчетному полю и, угловые

х

Полученные значения С представлены на рис. ?? (кривая 1). В проточной части водохранилища х < 103 км величина С составляет около 20. С приближением к озерной части скорость течения падает, а коэффициент Шези увеличивается до С ~ 30, а в наиболее широкой и глубокой части водоема числитель в расчетном выражении переходит через нуль и становится знакопеременным. Это означает, что в области 105 < х < 115 км решения уравнения не существует и понятие гидравлического трения, введенное в уравнениях Сен-Венана, становится условным. Причина этого феномена, возможно, в том, что уравнения мелкой воды, лежащие в основе системы Сен-Венана, выведены в приближении длинных волн^ когда глубина потока считается малой по отношению к горизонтальным масштабам волн. Резкое изменение морфометрических параметров при переходе к широкому плесу с увеличением глубин обусловливают возможность развития не только длинных, но и сравнительно коротких поперечных волн, что и приводит к

Рис. 2. Изменение коэффициента Шези С (кривая 1), параметра шероховатости п (кривая 2), средней скорости (кривая 3) и размеров поперечника (кривая 4) вдоль русла Новосибирского водохранилища

обнаруженному противоречию (относительное изменение ширины поперечников видно по кривой 4 на рис. ??).

Параметр шероховатости п в данном расчете на участке проточности составил величину ~ 0.03, которая вблизи критического участка, где С неограниченно возрастает, уменьшается до нуля. Значение п = 0.03 является характерным для водохранилищ равнинного типа.

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

2. Модель переноса взвеси

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

с 3 ? Сп частная

N

мутность отдельной фракции, причем считается, что с = ^ Сп, где N — число рас-

п=1

сматриваемых компонентов полидисперсной взвеси. Уравнение переноса концентрации получим из общего уравнения турбулентной диффузии твердой примеси [?] путем его осреднения по поперечнику

дЬс д Ьис д Ьтс дЬс д д с д д с

"ттт + ^----= т— Ькх— + — ЬКХ —, (10)

дЬ дх дх дх дх дх дх дх

Обозначим через Рх = Ь ( (т — тд)с — ) поток мутности в направлении х и

где тд — гидравлическая крупность наносов, вычисляемая по формуле Стокса [?].

дс ' дх/

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

Рх = 0 при х = (. (11)

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

^ = <>• (12)

дх

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

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

Второе обстоятельство связано с физической неполнотой условия (??), которое описывает лишь процессы гравитационного оседания (Рг < 0) и не учитывает процессы взмучивания мелких фракций, т. с. вовлечения наносов со дна в турбулентное движение. Это неоправданно с точки зрения правильного описания массообмена между дном и потоком, так как в природных водоемах осаждение в чистом виде наблюдается лишь при малых скоростях течения.

В неравномерных потоках с переменными скоростями течения приемлемым является условие

с = сь при х = (13)

где сь — значение мутности в данной точке дна. Для расчета величины сь используются фактическая информация о средней мутности потока на входном створе за исследуемый период и теоретическая формула Г. Эйнштейна [?] для расхода донных наносов п~и фракции

дш ~ Рп, (14)

р 1 п

I А ОПР Рш З^П I

ехр ( 0.39--— ) — 1

рт п*

п* — динамическая скорость, (п = 1,...,М) — диаметры частиц, р3 — плотность взвеси. Следуя [?], предположим, что распределение донной концентрации подчиняется зависимости, подобной (??), и запишем

СЬп = СоРп, (15)

где с0 — нормирующая величина, определяемая из условия равенства средней расчетной мутности фактическому значению

N

^ ^ Рпсп с£ай ,

п=1

с

С

„ - 1 [ „

данных наблюдении, сп — ~г I сп(х) средняя расчетная концентрация, найден-

ь ]

ная из решения задачи с краевым условием (??)■ Профиль сп(х) служит граничным значением концентрации при х = х\, а для его расчета используется одномерный стационарный аналог уравнения (??)■ На выходном сечении задавалось

д с

т— = 0 при х = х2. (16)

дх

Данные о мутности р. Обь [?] показывают, что в разные годы и от сезона к сезону концентрация взвешенных наносов в районе г. Камень-на-Оби (начальный створ Новосибирского водохранилища) испытывает значительные вариации, изменяясь в пределах 200...1000 г/м3, причем минимальные значения соответствуют осенне-зимней межени, а максимум твердого расхода наолюдается в весеннее половодье. Рассматривая режим мутности при расходе воды Q1 = 3000 м3/с, зададимся средним значением концентрации = 400 г/м3. Анализ данных наблюдений показал, что взвесь в районе г. Кам-ня-на-Оби содержит частицы диаметром 0.001...0.1 мм, причем качественно процентное соотношение фракций остается примерно стабильным в течение сезона. Для проведения конкретного расчета взят гранулометрический состав, отвечающий весеннему стоку 1993 г. с тттестикомпонентной градацией по фракциям.

Рисунок 11. а иллюстрирует пространственное распределение вертикального потока взвеси Рг, полученное в результате интегрирования уравнения (??) на установление. Темные участки на рисунке соответствуют Рх > 0, т. е. переносу взвеси в верхние слои, светлые — Рх < 0, т. е. процессам осаждения (осветление воды), нулевой фон с Рг ~ 0 помечен литерой "0".

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

Рис. 3. Расчетное поле вертикального потока взвеси Рг (светлые участки соответствуют осаждению, темные — взмыву) (а) и горизонтальное распределение придонных значений Рх (мм/год) для расходов воды $1=3000 (кривая 1) и $1=2000 м3/с (кривая 2) (б). Кривая 3 показывает суммарное по сечению распределение взвеси

со дна. При этом формируется равновесная концентрация, отвечающая транспортирующей способности потока [?]. Ниже по течению в области 5 < х < 15 км вследствие торможения уменьшается и транспортирующая способность, что обусловливает формирование зоны осветления. Далее следуют области взвешивания, которые перемежаются с областями аккумуляции наносов. В озерной части водоема при х > 80 км доминирует слабая седиментация, и лишь вблизи выходного створа за счет сужения русла и ускорения потока развиваются процессы взвешивания наносов.

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

~ Рг

показан на рис. ??, б (кривая 1) в терминах удельной величины Рг = —-, переведенной

РзЬ

для наглядности в единицы мм год. Процессы осаждения наиболее интенсивны в рус~ ловой части водохранилища с максимумом абсолютной величины около 4 мм в год (при х = 35 км). Придонный поток массы Рг для меньшего расхода воды показывает кривая 2 на рис. ??, б. Области экстремальных значений практически не меняют своего положения при уменьшении расхода, а области слабоинтенсивных процессов при х = 45 км смещаются против потока. Кривая 3, отражающая суммарное содержание взвешенного аллювия в воде (Ьс), позволяет судить о степени осветления воды при транзитном переносе: если не учитывать слабый взмыв вблизи створа ГЭС, то потеря массы наносов в котловине верхнего бьефа за счет осаждения составляет около 85 %. Хотя эта величина и находится в согласии с натурными данными, полученная скорость годовой аккумуляции наносов является оценочной. Для расчета фактических объемов заиления требуется учет внутригодовых вариаций расходов воды и гидромеханических параметров наносов за репрезентативный период или за все время существования водохранилища.

Заключение

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

Автор благодарит A.A. Атавина за полезные обсуждения постановочной части работы.

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

[1] Барышников Н.Б., Попов И.В. Динамика русловых потоков и русловые процессы. Л.:

Гидрометеоиздат, 1988. 455 с.

[2] Salomon J.С. Modelling turbidity maximum in the Seine estuary // Ecohydrodynamic Proc. 12-th Int. Liege Colloq. Ocean Hydrodyn. Amsterdam, 1981. P. 285-317.

[3] Васильев О.Ф., Воеводин Л.Ф.. Никифоровская B.C. Численное моделирование стратифицированных течений в системах открытых русел и водоемах разветвленной формы // Вычисл. технологии. 2004. Т. 6, № 1. С. 26-41.

[4] Роди В. Модели турбулентности окружающей среды // Методы расчета турбулентных течений / Под ред. В. Кольмана. М.: Мир, 1984. С. 227-322.

[5] Шлычков В.А. Численное моделирование речных потоков с учетом генерации вихрей на границе русло—пойма // Водные ресурсы. 2008. Т. 35, № 5. С. 546-553.

[6] Квон В.П., Квон Д.В. Математическое моделирование процесса заглубления перемешанного слоя в стратифицированной жидкости // Прикл. механика и техн. физика. 1997. Т. 38, № 1. С. 82-85.

[7] Россинский К.П., Дебольский В.К. Речные наносы. М : Наука, 1980. 214 с.

[8] Гришанин К.В. Основы динамики русловых потоков. М.: Транспорт, 1990. 320 с.

[9] Белолипецкий В.М., Генова С.Н. Вычислительный алгоритм для определения динамики взвешенных и донных наносов в речном русле // Вычисл. технологии. 2004. Т. 9, № 2. С. 9-25.

[10] Зубкова K.M. Формирование и балане взвешенных наносов р. Оби на участке от г. Барнаула до г. Камень-на-Оби // Тр. Гос. гидрологического ин-та. 1991. № 349. С. 31-45.

[11] Кантаржи И.Г., Анцыферов С.М. Моделирование взвешенных наносов под волнами на течении // Океанология. 2005. Т. 45, № 2. С. 173-181.

Поступила в редакцию 13 мая 2009 г., с доработки — 18 июня 2009 г.