ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKII REGION. NATURAL SCIENCE. 2019. No. 3
ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ _PHYSICAL AND MATHEMATICAL SCIENCES_
УДК 532.5, 556 DOI 10.23683/0321-3005-2019-3-4-14
МОДЕЛЬ СТАЦИОНАРНОГО ТУРБУЛЕНТНОГО ТЕЧЕНИЯ И ПРОЦЕСС СЕДИМЕНТАЦИИ ПРИМЕСИ*
© 2019 г. М.Ю. Жуков12, Е.В. Ширяева1, А.В. Васильев1
1Южный федеральный университет, Ростов-на-Дону, Россия, 2Южный математический институт, Владикавказ, Россия
THE MODEL OF STATIONARY TURBULENT FLOW AND SEDIMENTATION PROCESS
M.Yu. Zhukov12, E.V. Shiryaeva1, A.V. Vasiliev1
1Southern Federal University, Rostov-on-Don, Russia, 2Southern Institute of Mathematics, Vladikavkaz, Russia
Жуков Михаил Юрьевич - доктор физико-математических наук, профессор, заведующий кафедрой вычислительной математики и математической физики, Институт математики, механики и компьютерных наук им. И.И. Воровича, Южный федеральный университет, ул. Мильчакова, 8а, г. Ростов-на-Дону, 344090, Россия; профессор, Южный математический институт, ул. Маркуса, 22, Владикавказ, РСО-А, 362027, Россия, e-mail: [email protected]
Ширяева Елена Владимировна - кандидат физико-математических наук, доцент, кафедра вычислительной математики и математической физики, Институт математики, механики и компьютерных наук им. И.И. Воровича, Южный федеральный университет, ул. Мильчакова, 8а, г. Ростов-на-Дону, 344090, Россия, e-mail: [email protected]
Васильев Андрей Владимирович - магистрант, кафедра вычислительной математики и математической физики, Институт математики, механики и компьютерных наук им. И.И. Воровича, Южный федеральный университет, ул. Миль-чакова, 8а, г. Ростов-на-Дону, 344090, Россия, e-mail: avvasi-lievl995@gmail. com
Michael Yu. Zhukov - Doctor of Physics and Mathematics, Professor, Head of Department of Mathematics and Mathematical Physics, Vorovich Institute of Mathematics, Southern Federal University, Milchakova St., 8a, Rostov-on-Don, 344090, Russia; Professor, Southern Institute of Mathematics, Markusa St., 22, Vladikavkaz, RNO-A, 362027, Russia, e-mail: [email protected]
Elena V. Shiryaeva - Candidate of Physics and Mathematics, Associate Professor, Department of Mathematics and Mathematical Physics, Vorovich Institute of Mathematics, Southern Federal University, Milchakova St., 8a, Rostov-on-Don, 344090, Russia, e-mail: [email protected]
Andrew V. Vasiliev - Magistrant, Department of Mathematics and Mathematical Physics, Vorovich Institute of Mathematics, Southern Federal University, Milchakova St., 8a, Rostov-on-Don, 344090, Russia, e-mail: avvasi-liev1995@gmail. com
На основе общих уравнений движения жидкости построена асимптотическая модель двухслойного пространственно одномерного стационарного турбулентного течения вязкой несжимаемой жидкости в области, бесконечной в горизонтальном направлении. В качестве малого параметра при построении главных членов асимптотических разложений используется отношение характерных размеров в вертикальном и горизонтальном направлениях. Кинематическая турбулентная вязкость задается кусочно-гладкой, постоянной в верхнем слое и зависящей от координат в
* Работа выполнена при финансовой поддержке базовой части государственного задания № 1.5169.2017/БЧ Министерства образования и науки РФ, ЮФУ.
ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKII REGION. NATURAL SCIENCE. 2019. No. 3
нижнем. Кроме этого, считается, что на верхней границе вязкость обращается в нуль. Особенностью модели является отказ от условий прилипания вязкой жидкости к нижней границе нижнего слоя. Показано, что в окрестности нижней границы, в тех её частях, которые соответствуют впадинам (ямам), имеются придонные вихри. В указанном стационарном течении исследуется квазистационарный процесс седиментации пассивной взвешенной примеси. Профили течения жидкости определяются аналитически. Для численного решения задачи о седиментации примеси в поле тяжести использован метод конечных элементов, позволяющий эффективно решать задачи с кусочно-гладкими коэффициентами. Структура построенного аналитического решения для течения жидкости позволяет при дополнительных предположениях без труда использовать и классический вариант краевых условий, соответствующих прилипанию вязкой жидкости к нижней границе или её фрагментам.
Ключевые слова: вязкая жидкость, турбулентная вязкость, седиментация, придонные вихри, метод конечных элементов.
On the basis of general equations offluid motion the asymptotic model of two-layer spatially one-dimensional stationary turbulent flow of viscous incompressible fluid in the region infinite in the horizontal direction is constructed. The coefficient wich is characterized scales in the vertical and horizontal directions is used as a small parameter in the construction of the principal terms of asymptotic expansions. The kinematic turbulent viscosity is defined as piecewise smooth-constant in the top layer and dependent on the coordinates in the bottom layer. In addition, we assume that viscosity is zero at the outer boundary of the top layer. A feature of the model is the rejection of the non-slip boundary condition at the lower boundary of the bottom layer. It is shown that in the vicinity of the lower boundary there are near-wall vortices in those parts of boundary that correspond to the non-straight fragments with negative curvature. The quasi-stationary process of the passive impurity suspended sedimentation is investigated in this stationary flow. The fluid flow profiles are determined analytically, and for the numerical solution of the sedimentation problem in the gravity, the finite element method is used, that allow us effectively to solve the problems with piecewise smooth coefficients. The structure of the analytical constructed solution for the fluidflow allows us to easily use the classical version of the non-slip boundary conditions at the lower boundary or its fragments.
Keywords: viscous liquid, turbulent viscosity, sedimentation, bottom vortices, finite element method.
На основе асимптотической модели двухслойного пространственно одномерного стационарного турбулентного течения исследуются седиментация и перенос пассивной примеси в бесконечных горизонтальных слоях. Разделение течения на два слоя является по сути фиктивным и связано с заданием коэффициента кинематической турбулентной вязкости в виде непрерывной кусочно-гладкой функции вертикальной координаты с разрывом производной на границе между слоями. Такой выбор вязкости, в частности, предложен в монографии [1, с. 33-35] и достаточно хорошо согласуется с известными экспериментальными измерениями в реальных течениях жидкости [2, гл. 20]. В нижнем слое турбулентная вязкость выбрана зависящей от вертикальной координаты. Горизонтальная компонента скорости имеет логарифмический профиль, типичный для турбулентных течений. В верхнем слое вязкость не зависит от вертикальной координаты. Горизонтальная скорость имеет степенной профиль, характерный для ламинарных течений. На границе между слоями вязкость непрерывна. Описанный выше вид вязкости предложен в [1] для однослойной модели с плоскими границами горизонтального слоя — верхней свободной и нижней плоской (дно водоема). В середине слоя вязкость непрерывная, гладкая и имеет разрыв второй производной. Для предлагаемой модели приведены различные гидрологические соотношения и обоснования пре-
имуществ такого выбора. Представляется более естественным описывать течение с указанной вязкостью при помощи двухслойной модели — в каждом вертикальном сечении в верхнем слое вязкость постоянная, в нижнем зависит от вертикальной координаты. Анализ показал, что модель двухслойного течения, в частности, избавлена от основного недостатка моделей для однослойного течения — необходимости использования гипотезы о «мгновенной адаптации», по которой профиль свободной поверхности жидкости полностью повторяет рельеф дна водоема (по поводу критики такой гипотезы см. [3-7]).
Асимптотическая модель, являющаяся вариантом гидростатической модели для стационарного течения однослойной жидкости с постоянной вязкостью, предложена в работе [8], где обоснована гипотеза о «мгновенной адаптации» и приведено численное исследование поведения пассивной примеси конечно-разностными методами. Предлагаемая ниже асимптотическая модель двухслойного течения принципиально не отличается от модели в [8] в том смысле, что использованы те же самые предположения о порядках малости параметров задачи. Однако выбор кусочно-гладкой турбулентной вязкости приводит к решению, которое существенно отличается от указанного в [8].
Особенностью предлагаемой модели двухслойного течения является замена классических крае-
ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKII REGION.
NATURAL SCIENCE.
2019. No. 3
вых условии прилипания на нижнеи границе кинематическим условием, характерным для свободной деформируемой поверхности. Для описания стационарного течения жидкости с профилем горизонтальной скорости, логарифмическим для нижнего слоя и степенным для верхнего, построено аналитическое решение, анализ которого позволил для некоторых значений параметров обнаружить наличие придонных вихревых структур. Оказалось, что такая модель дает возможность задавать произвольные профили границ — нижней, между слоями и свободной верхней. Структура предлагаемого решения при некоторых предположениях позволяет легко перейти к решению для классических условий прилипания.
Задача о седиментации и переносе пассивной примеси решается численно при помощи метода конечных элементов. Кинетические коэффициенты уравнения диффузии с переносом, в отличие от модели, предложенной в [8], выбраны более реалистичными. Коэффициент диффузии не считается постоянным. Он пропорционален вязкости (в соответствии с так называемой аналогией Рейнольдса [9]). Для скорости седиментации примеси в поле тяжести выбрана простейшая формула Стокса, соответствующая шарообразным частицам. Ее без труда можно заменить любой из предлагаемых в многочисленных работах формул [3, с. 51, 154-158], учитывающих иные формы частиц, их размеры и т. д. Выбор метода конечных элементов для решения обусловлен, в частности, тем, что в его основе лежит слабая формулировка задачи, позволяющая без труда строить решения задач с кусочно-разрывными коэффициентами (в данном случае турбулентной вязкости и диффузии). Это дает возможность решить задачу о переносе примеси для всей области в целом, а не для каждого слоя области в отдельности. Специальные приемы аппроксимации уравнений позволяют существенно повысить точность выполнения уравнения баланса массы.
Основные уравнения и соотношения
Рассматривается седиментация примеси в стационарном пространственно двухмерном турбулентном течении несжимаемой вязкой жидкости в случае, когда кинематическая турбулентная вязкость задается соотношениями:
= (1)
цг(х) = (п(х) — ЛЬ(х))(пЧх) - л(х)),
^ь(х,г) = (г — (х))(т^г (х) — г).
Здесь Кг - характерная безразмерная турбулентная вязкость; ц1(х), ць(х), ц(х) - функции, задаю-
щие свободную границу жидкости, рельеф дна и некоторую промежуточную границу, на которой происходит изменение характера поведения вертикального профиля вязкости (в части формул (1) аргумент х опущен для сокращения записи); х, z - горизонтальные и вертикальные координаты.
Соотношения (1) для случая плоских границ (rj1,
h
r¡u, ц - постоянные, зависимость вязкости от горизонтальной координаты х отсутствует) предложены в [1, с. 33-35] и достаточно хорошо коррелируют с результатами экспериментов [2, гл. 20; 10].
Наличие поверхности z = tf(x), на которой профиль вязкости по переменной z непрерывен (^(х) = /ль (x,ij(x))) и в общем случае кусочно-гладкий (^(х) ф (х, ij(x))), делает разумным рассмотрение модели двухслойного течения жидкости.
Для описания поведения жидкости и примеси используем упрощенные асимптотические уравнения, записанные в безразмерном виде в прямоугольной системе координат. Ось х наклонена к горизонтали под углом а, направление оси z (при а = 0) противоположно направлению действия силы тяжести: Ux + wz = 0, (2)
tz + g0 sin а = 0, т = p0uz, (3)
0 = -pz- до cos а, (4) ct + исх + wcz + iz = 0, (5)
1 = —D0cz + cS0 cos a. (6) В (2)-(6) u(x, z), w(x, z) - горизонтальная и вертикальная компоненты скорости; т(х, z) - касательное напряжение; р(х, z) - давление; с(х, z, t) - концентрация примеси; i - нормальная компонента плотности потока примеси; D0 (х, z) - коэффициент диффузии; S0 (х, z) - скорость седиментации (S0 < 0 соответствует движению примеси в поле тяжести по направлению ко дну водоема); д0 - величина ускорения силы тяжести; (х, z) - кинематическая вязкость жидкости; а - угол наклона к горизонтальной плоскости.
Уравнения справедливы для верхнего = = {(х,г)\^1(х) < z < ц(х)} и нижнего = = {(x,z): ij(x) < z < r¡b(x) + z0(x)} слоев, где z0(x) - шероховатость границы z = tfb(x), которую необходимо рассматривать в теории для того, чтобы при выбранной вязкости (1) не возникало сингулярности скорости течения на границе (рис. 1) [2, гл. 20]. В дальнейшем величины, относящиеся к различным слоям, снабжены верхними индексами: для верхнего слоя (t) и для нижнего (Ь).
Простейшие асимптотические уравнения, аналогичные (2)-(6), рассматривались и численно исследовались для случая одного слоя и постоянных коэффициентов вязкости и диффузии в [8], более общие двухслойные модели, описывающие нестационарные процессы седиментации примеси, - в [3-7].
ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKII REGION.
NATURAL SCIENCE.
2019. No. 3
Заметим, что (2)—(6) являются упрощением уравнений классической гидростатической модели - давление определяется уравнением гидростатики (4), но, кроме этого, при определении горизонтальной скорости с помощью уравнения (3) пренебрегают инерциальными членами и градиентом давления в горизонтальном направлении.
Асимптотические уравнения получены на основе общих уравнений, описывающих поведение вязкой несжимаемой жидкости, движущейся по наклонной плоскости под действием силы тяжести и содержащей пассивную примесь. В безразмерных переменных они имеют вид
div V = 0, = д + V ■ V,
dt dt
dv
P~¿ = —Vp + div T + pf, T = p^(Vv + (Vv)T),
de
P~t + div i = °, (7)
i = —pDVc + pcS, S = S(— sin a, cos a), F = g (sin a, —cosa).
Здесь v = (и, w), и, w - горизонтальная и вертикальная компоненты скорости; р - давление; Т -тензор вязких напряжений; с - концентрация примеси; i - плотность потока примеси; D - коэффициент диффузии; S - скорость седиментации (S < 0 соответствует затоплению примеси в поле тяжести); р - плотность жидкости; д - величина ускорения силы тяжести; ^ - кинематическая вязкость жидкости; а - угол наклона плоскости к горизонтальной поверхности.
Связь между безразмерными и размерными переменными, отмеченными тильдой, дается соотношениями:
(t,x,ü) = (T,t,L,x,L,T, 1u), (с,р,ß,§У=(с,с,р,р,^v, д,д), (г,ц1,ц,ць,г0) = H„(z,rf ,-ц,-ць ,z0), (w,S) = H,T,-1(w,S0), р = p,L2,T~2p, р = 1,
= ßT,L,2, Dg£2 = DTtL-
9o = zgtT^Lj,1,
1-2
(8)
(9)
- äi = L,'
частности, уравнения (2)-(6), справедливы для достаточно широкой области изменения параметров, необязательно совпадающей с той областью, которая требуется для строгого математического вывода.
Гидрологические определяющие соотношения.
Система уравнений (2)-(6) не является замкнутой в том смысле, что в отличие от коэффициента турбулентной вязкости ^0 (1) конкретный вид коэффициента диффузии О0 и скорости седиментации 50 не указан.
Коэффициент турбулентной вязкости. Для размерной величины коэффициента турбулентной вязкости в случае установившегося турбулентного течения при исследовании русловых потоков обычно используется соотношение [2, гл. 20; 11, (4.10); 12, (42)]
Д = ку.Н.^ь(х,г), (11)
где к - постоянная Кармана (« 0,4); V. - динамическая скорость течения в русловом потоке; Н. - глубина потока; - безразмерная функция, определенная формулой (1).
В случае использования соотношения (11) размерная горизонтальная скорость течения й определяется формулой [11, (4.2); 12 (42.4)]
и=~Ы~ (12)
При определении динамической скорости V. в случае открытого русла применяется формула [11, (3.9)]
2
g,H,I, I & а & sin а,
(13)
где I - гидравлический уклон.
Для согласования (11)—(13) с использованным переходом к безразмерным величинам (8), (9) следует выбирать величину динамической скорости V. в качестве характерной скорости и параметр Кг, входящий в (1), в виде
Kz = к, v,=
(14)
Обратим внимание на то, что уравнение (3) после подстановки формулы (1) содержит безразмерный параметр
до sin а 1 д„Н„ sin а
Здесь Т., р., е., д. - характерные величины времени, плотности, концентрации, ускорения силы тяжести; Ь., Н. - характерные размеры в направлениях х, г; £ - отношение характерных размеров.
Уравнения (2)-(6) являются главным членом асимптотического разложения уравнений (7) при £ ^ 0 и при выполнении условий
до = 0(1), ро = 0(1), (10)
0о = 0(1), Бо = 0(1).
Соблюдение соотношений (1 0) не является критическим. Как правило, асимптотические модели, в
Kz к KV2
(15)
который в [8] обозначается 2Я и трактуется как число Рейнольдса. Действительно, в случае ламинарных течений величина Кг аналогична обычной вязкости (не определяется формулой (1 4)), и такая трактовка разумна. Однако в случае турбулентных течений с универсальным логарифмическим профилем указанный параметр также является универсальным и равным к \
Коэффициент турбулентной диффузии. Для коэффициента турбулентной диффузии используем
ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKII REGION.
NATURAL SCIENCE.
2019. No. 3
аналогию Рейнольдса - пропорциональность коэффициента турбулентной диффузии и кинематической турбулентной вязкости [2, с. 269; 9, с. 287].
О0(х,г) = 8^0(х,г), (16)
где 8-1 = Рг - диффузионное число Прандтля (число Пекле Ре = (ек5)-1).
В гидрологии под аналогией Рейнольдса иногда понимают равенство всех размерных кинетических коэффициентов - турбулентной кинематической вязкости, турбулентной теплопроводности и, в частности, турбулентной диффузии (¡1 = 0 [11, (3.29)]). Однако, как указано в [11, с. 48], для турбулентной диффузии такая аналогия работает лишь при малых значениях концентраций и для мелких взвешенных частиц примеси. Именно по этой причине в формуле (16) присутствует дополнительный параметр - диффузионное число Прандтля.
Скорость седиментации. Для скорости седиментации частиц примеси используем формулу Стокса, определяющую скорость движения сферических частиц в вязкой жидкости [11, 12], которую удобно записать, используя размерные величины:
(17)
где - характерный диаметр частиц примеси; р.* -обычная (нетурбулентная) размерная кинематиче-
с
ская вязкость жидкости; р; - плотность частиц; р* -плотность жидкости.
В речных потоках плотность твердых частиц р1 изменяется в достаточно узких пределах 26002700 кг м-3. Характерный размер частиц й* в зависимости от их типа можно, например, задавать при помощи таблицы из [11, с. 186; 3, с. 147].
Формулу (17) с учетом (14) удобно представить в форме
50 (х, z)cosa = в0, (18)
Sn =
g,d.ÍL, cosa ,
osa / _ p£\
v, V pj'
18ц*Н*у* \ р*/ Взамен (1 7) можно использовать любую из многочисленных формул, имеющихся в литературе (монография [3, с. 51, 154-158] и цитируемые в ней источники). Предписываемый формулами (10) порядок величины 50 = 0(1) легко обеспечить выбором величин й*, р*. Это можно трактовать и как некоторое ограничение на размер и плотность частиц, для которых применимо асимптотическое приближение. Оно, как вновь уместно отметить, в действительности может иметь более широкую область применимости, чем диктуемую ограничениями (1 0). Обратим внимание на то, что в формуле (18) присутствует обычная вязкость р*, а не турбулентная р.. Иными словами, турбулентный механизм процесса переноса примеси в вязкой жидкости под действием силы тяжести не предполагается.
Шероховатость. Параметр шероховатости z°(x) (точнее, функция) является феноменологическим параметром. Существует большое количество теорий и полуэмпирических соотношений для его определения [1, формулы (2.25)-(2.27); 3, с. 147— 150; 11, с. 39, 40]. Хотя теория и позволяет считать шероховатость зависящей от координаты х, чаще всего выбирается z°(x) — const. Далее считаем [3, с. 150, табл. K.5; 13]
— 15 , z° — 0,033 kg, где k* - эффективный выступ шероховатости; d% -характерный размер шероховатости.
Модель двухслойного течения жидкости
Уравнения (2)-(6), описывающие процесс седиментации примеси, расщепляются - предварительно из (2), (3) определяется профиль стационарного течения жидкости, затем из (5), (6) — нестационарный профиль концентрации.
Уравнения (2), (3) следует дополнить краевыми условиями. Удобнее решать задачу и формулировать краевые условия при помощи функции тока, позволяющей автоматически удовлетворить уравнение (2). Для каждого слоя функция тока задается соотношениями:
и* — té, w* — (19)
ub — i¡}b, wb — -грь.
Полная постановка задачи для определения профилей стационарного течения с учетом (19) имеет вид (у функций rf(x), r](x), r]b(x) для сокращения записи опущены аргументы)
t¡ + g0 sin а — 0, т1 — Кгр1-ф1гг, rb+gosina — 0, ть — Kzpb^bz, (20)
Шх,Л*) — 0, (rt(x,r1t) — 0), (21)
^t(x,r¡t) — Qt + Qb, (22)
^t(x,r]) — Qb, фb(x,r]) — Qb, (23)
/(*Wzz(x,r¡) — pb (x, r¡)ipbz(x, r¡), (24)
фb(x,r]b(x) + z0(x)) = 0,
(25)
V1
Qz = Iv uz(x,z) dz =ipt(x,rt) -^z(x,r),
Qb = tfb+zoUb(x,z)dz= (26)
= \pb(x,r)- ^b (x, r\b + z°).
Здесь Qb, Ql - заданный расход жидкости в каждом слое (Qb, Ql - const). Динамические условия (21), (24) следуют из непрерывности на границах нормальной компоненты полного тензора напряжений ([-рп + Т • п] = 0, где [•] означает величину разрыва функции) и соответствуют отсутствию касательного напряжения т на свободной границе и непрерывности касательного напряжения между слоями. Кинематические условия (22), (23), (25) эквивалентны отсутствию на границе нормальной
ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKII REGION.
NATURAL SCIENCE.
2019. No. 3
компоненты скорости (V • п = 0). Два кинематических условия (23) возникают ввиду неподвижности границы 2 = ](х) между слоями И,', И,й.
Интегрируя уравнения (20) по г и используя динамические условия (21) и (24), получим с учётом (15)
к
tit t nt(x)-z
(27)
+
6к^(х)
(г-у{ХУ-У{)2 + z-n Qt + Qb
6Kßt(x)
Tpb(x,z) =
(z-nЬ)Ы*-УЬ
(28)
(29)
(n-nb)(z-nb-z°)]nn-nb ,
+ Q b
В асимптотической модели такое условие эквивалентно отсутствию касательной компоненты скорости. Условие (25) для функции тока сохраняется и соответствует непроницаемости границы. Это означает, что функция тока может быть получена из соотношения (29) при дополнительном условии (30) и имеет вид
Tpb(x,z) =
ln
z-n"
z-nb-z0
Обратим внимание, что в силу краевых условий (21), (24) соотношения для касательных напряжений ть, т1 в слоях И,й, И,' совпадают между собой.
Решение задачи (27), (22), (23), (25)-(26) записывается в форме
Такое решение существует лишь в случае, когда ~ц(х) — -ць(х) = h = const, (31)
Qb = const, Qz = const, z0(x) = const, где h - постоянная толщина слоя Lb.
Функция тока по-прежнему определяется формулой (28), величину Qz можно задавать произвольной константой, а величина Q b уже не является произвольной константой и определяется соотношением
" ' ' * (32)
к (т]-т]ь-г0) (г}-г]ъ-г°)-
Величины QQb можно задавать произвольными константами. Скорости и', шиь, шь легко вычисляются при помощи (19) и не приведены ввиду громоздкости.
Формально решения (28), (29) справедливы при любых достаточно гладких функциях ](х), ](х), ]ь(х), г0(х). Однако на практике естественно требовать, чтобы при |х| ^ от указанные функции стремились к некоторым константам, т. е. горизонтальные слои !,й, являлись плоскими на бесконечности. Это автоматически приведет к стремлению на
бесконечности горизонтальных компонент скорости
Ь ' ь
и", и' к некоторым константам, а вертикальных ш ,
- к нулю.
Анализ решения. Кинематическое краевое условие (25), соответствующее непроницаемости для жидкости границы г = ]ь(х) + г0(х) (т.е. V • п = 0), не совпадает с обычно используемым условием прилипания (V = 0). Однако граница г = ]ь(х) + г0(х), строго говоря, не является твердой - это всего лишь (фиктивная) граница шероховатостей, введенная в рассмотрение с целью ликвидации сингулярностей скорости вблизи границы г = ] (х). Более того, при достаточно малых шероховатостях г0(х) ^ +0 вязкость жидкости ^о(х) стремится к нулю (1), и требование прилипания жидкости к границе при фактически исчезающей вязкости ничем не обосновано.
Выполнение условия прилипания влияет на характер стационарного течения. Действительно, пусть выполнено условие прилипания
иь(х,-пь+ г0) = -ф%(х,-пь+ г0) = 0. (30)
Qb = Qb=±.(hln±-(h-z0)).
Здесь обозначение Q0) введено для удобства.
Соотношение (31) показывает, что для нижнего слоя в случае выполнения условия прилипания (30) имеет место эффект мгновенной адаптации -профиль границы х = ](х) отличается от профиля границы 2 = ]ь(х) на постоянную величину. Мгновенная адаптация типична для однослойных моделей, в которых свободная верхняя граница повторяет рельеф дна. Преимуществом двухслойной модели является тот факт, что верхняя свободная граница слоя И,' может быть произвольной и не обязана повторять рельеф дна, что, на наш взгляд, является более приемлемым. Следует отметить, что и соотношение (32) указывает на недостаток модели с прилипанием - при г0(х) ^ +0 расход Qb ^ от, что, конечно же, является следствием возникновения синь
гулярности для скорости и".
В пользу отказа от условия прилипания (30) и сохранения лишь условия непроницаемости границы (25) говорит и тот факт, что все функции ] (х), т](х), ]ь(х), определяющие границы, могут быть заданы произвольно. Имеется возможность моделирования любых ситуаций, возникающих на практике. Строго говоря, любые ограничения, в частности (31), следует считать «дефектом» модели.
Обратим также внимание на то, что выбор условия прилипания никак не влияет на величину касательного напряжения т = к-1(т]' —г]ь — г0) (27) на нижней границе 2 = ]ь(х) + г0(х), и отказ отусло-вия прилипания можно трактовать как выполнение на границе условия скольжения Навье - касательная скорость пропорциональна касательному напряжению. В современных работах условие Навье объясняется при помощи введения слоя проскальзывания
0
к
z
к
к
ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKII REGION.
NATURAL SCIENCE.
2019. No. 3
[14; 15; 16, с. 37-40], который для рассматриваемой модели соответствует слою с шероховатостью.
Наличие в решении (28), (29) двух независимых констант ((1 (и особенно то, что константа на самом деле в случае прилипания определяется соотношением (32)), конечно, выглядит неестественно. Требуется дополнительное соотношение, позволяющее связать эти константы. Наиболее разумным является выбор совпадения средних горизонтальных скоростей течения жидкости в нижнем и верхнем
слоях ¿ау(
(33)
Здесь и1У, и0у - средние скорости, которые опре-
деляются соотношениями
Qz
ulv(x) =
rjt(x)-rj(xу Qb
(34)
ub = —
av r(x)-rb(x)-z°(xy
где поток концентраций (6) и скорости и, ш определены с учетом соотношений (1), (14), (16), (18), (19), (28), (29)
8^о(х,г)
=
( u
■cz + S 0С,
(uf, wf), (х, z) е Ii,
Ввиду того, что и0у, щу зависят от х, равенство средних скоростей возможно, лишь когда границы являются плоскими, что и отражено соотношением (33), взамен которого можно использовать условия при х ^ +от, но естественнее задавать такие условия на входе потока жидкости, т.е. при х ^ —от.
Обратим внимание на то, что для модели (20)-(26) (т.е. без требования условия прилипания) при достаточно больших величинах полного расхода жидкости ( = ( + возможно возникновение противотока, когда в некоторой части области горизонтальные компоненты скорости иь (х, г) будут отрицательными [17, с. 11]. В частности, нетрудно получить условие, при выполнении которого горизонтальная скорость на границе и0 (х,ць (х) + г0(х)) будет неотрицательна. Если использовать (33), (34), то
пЬ _ Уо-Уо-го п ( = Vto-Vbo-z0(, _ _По—По_п ( = Vtt-Vbo-z0(,
где ц0, ц0, г]0, - значения функций при х ^ —от.
Неравенство, соответствующее неотрицательности горизонтальной скорости (иь(х,ць(х) + +г0(х)) > 0) на нижней границе слоя выполнено лишь при условии (}ь > (д, где ((0 определено соотношением (32) при 1г = ц0 — ц0.
Движение примеси в стационарном потоке
Движение примеси в стационарном потоке, задаваемом (28), (29), в области и описывается уравнениями (5),
Сь + исх + шсг + = 0, (35)
(36)
0 шь), (х,г)£Ль. (37)
Систему (35)-(37) следует дополнить краевыми и начальными условиями. На верхней свободной границе слоя И,' выбираем естественное условие отсутствия потока концентрации
1(х,^) = 0. (38)
Такое же условие можно задавать и на нижней границе слоя. Формально при |х| ^ от задание каких-либо условий не требуется, однако при численном решении бесконечную область и приходиться заменять ограниченным в горизонтальном направлении конечным фрагментом. В этом случае для левой границы можно выбрать условие, аналогичное (38), для правой - условие полной проницаемости, которое соответствует сохранению концентрации примеси при переносе через границу и имеет вид [18, с. 113, 114]
Сь + исх + шсг = 0, х = хг, (39)
где хг - координата правой границы конечной области.
Более детальная информация о выбранных граничных условиях для примеси приведена ниже при описании результатов вычислений.
Уравнения (35), (36) имеют разрывные коэффициенты. Действительно, на линии г = ц(х), т.е. границе между слоями, имеются разрывы скоростей и, ш и производной вязкости (10 по переменной х 5и(х) = и1(х,-ц) — иъ(х,ц) Ф 0, (40)
8ш(х) = шг(х,-ц) — шь(х,ц) Ф 0, р\(х) Ф$(х,ц).
По этой причине для численного решения (35)-(37) удобно применять метод конечных элементов, который использует переход к слабой формулировке задачи, не чувствительной к конечным разрывам коэффициентов.
Обратим внимание на то, что х = ц(х) является также линией уровня функции тока, на которой нормальная компонента полного потока примеси = = \с + \ непрерывна ввиду того, что V • п = 0 (см. кинематическое условие (23)), а нормальная компонента потока • п по причине непрерывности вязкости /л,:(х) = /ль(х, ц). Это, в частности, означает, что если вместо объединения областей и решать задачу для каждой области и отдельно, то перенос примеси через границу между областями осуществляется непрерывным образом по нормали к границе, а разрывы горизонтальной скорости на границе будут приводить лишь к искажению профиля концентрации в окрестности границы в горизонтальном направлении.
ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKII REGION.
NATURAL SCIENCE.
2019. No. 3
В отличие от прямых конечно-разностных методов при использовании метода конечных элементов, как правило, трудно добиться численного выполнения закона сохранения (баланса) массы, в первую очередь в связи с трудностями выполнения уравнения неразрывности div V = 0. Оказалось, что при специальном разбиении конвективного члена в уравнении (35)
V -Чс -Чс + ¿гу(>с)), (41)
удается существенно повысить точность решения.
Результаты вычислений
Для численного исследования использован метод конечных элементов и пакет FreeFem++ [18, 19]. Расчеты приведены при значениях параметров, характерных для равнинных рек: ширина реки W* = = 100 м; глубина реки Н* = 10 м; характерный горизонтальный размер Ь* = 150 м; ускорение силы тяжести д* = 9,81 м • с-2; плотность жидкости р* = = 1000 кг • м-3; плотность частиц примеси р* = = 2650 кг • м-3; уклон а = 0,1256 - 10-3 м • км-1; характерный размер частиц примеси й* = 5,74х х10-5 м; характерный размер шероховатости й* = = 0,3 м; объемный расход Q0 = 680 м3 • с-1 (рас-
п ^
ход при заданной ширине реки Q*=^' =
= 6,8 м2хс-1). Кроме этого, задана постоянная Кармана к = 0,4 и число Прандтля З-1 = 20. В этом случае динамическая скорость V* = 0,11 м • с-1; характерная средняя скорость и* = — = 0,68 м • с-1;
н*
обычная кинематическая вязкость жидкости ^* = = 10-6 м • с-2; характерное время Т* = 1350 с. Тогда безразмерные расход жидкости Q, скорость седиментации 50, шероховатость г°(х) и параметр £
будут
Q = 6,125607624, б0 = -0,4, (42)
г0(х) = 0,015, е = 0,067.
Форму слоев жидкости определяем соотношениями:
= - а^Шх - Х1)), 2 (43)
V = Ч0 + а2(х - х2)е-^2(х-х2)2,
(42), соответствует случаю, когда иь(х,ць(х) + +г0(х)) = 0 при |х| ^ от. Иными словами, всюду, где границы слоев плоские, на нижней границе области выполнено условие прилипания, которое нарушается лишь вблизи отклонения формы границ от плоской.
Величины расходов жидкости в слоях Qb, Qt при выбранных параметрах, указанных на рис. 1, Qb = = 2,612, Qt = 3,514.
На рис. 1 показаны также линии уровня нескольких функций тока (фг(х,цг) = Qb + Qг, ^(х,ц) = Qb, ^ь(х,ч) = Qb, *рь(х,ць) = 0), очевидно, совпадающие с границами слоев г = цг(х), х = ч(х), г = чЬ(х). Стрелки на рис. 1 соответствуют скоростям (и, ш + з0), с которыми двигались бы частицы примеси при отсутствии диффузии. Ввиду разрывов скоростей течения (8и Ф 0, Зш Ф 0) частицы в разных слоях в окрестности границы имеют разные скорости движения (двойные стрелки на линии ^(х,ч(х)) = Qb). Это позволяет понять механизм седиментации примеси в окрестности границы между слоями. Двигаясь в верхнем слое по линии тока трг(х,ц(х)) = Qb, частицы примеси смещаются за счет скорости седиментации и диффузии на линию тока ^(х,ч(х)) = Qb, продолжая движение в нижнем слое ИЛ
ГЬ = rb + а3(х - х3)е-^з(х-хз)2 Здесь функция гЬ(х) задает рельеф дна в форме = 0,25; х2 = 1,5; Р2 =2;Ло = -1; аз = 0,5; хз = 1;Рз = 2
«яма-горб» с размером а3, ^(х) определяет форму
свободной поверхности с перепадом высот а1, т](х) -
форму границы между слоями; параметры а^, р, х^
ь г
задают амплитуды, сдвиги и сглаживания; ч0, ч0, ч0, определяют расположения границ слоев при х ^ -от (значения параметров даны на рис. 1).
Обратим внимание на то, что для приведенных результатов расчетов величина ч0, задающая положение границы между слоями и профили границ
Рис. 1. Расположение слоев и функции тока ipb (х, z), фг(х,г). Функции rb, Г, Г заданы соотношениями (43); Г0 = -0,1; а1 = 0,1; х1 = 0; р1 = 1, r0 = -0,565; а2 =
= 2; rb = -
/ Fig. 1. Layers location and stream functions ipb(x,z), фг(х,г).
Functions rb ,Г,Г are given by relations (43); rO = -0.1; а1 = 0.1; x1 =0;P1 = 1; r0 = -0.565; а2 = 0.25; x2 = 1.5; P2 = 2;r0 = -1; а3=0.5;х3 = 1;Р3=2
Пунктир на z = гЬ(х) + г0(х) указывает области, для которых в окрестности линии горизонтальные скорости отрицательны (ub < 0). В таких областях либо возникает противоток жидкости [17,
ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKII REGION.
NATURAL SCIENCE.
2019. No. 3
с. 11], либо, как в данном конкретном случае, происходит образование придонных вихрей, которые показаны на рис. 2.
Для расчета структуры вихрей и границ областей, в которых они расположены, использовались уравнения
— = мй(х,г), — = шй(х,г)
^ 4 у Л: у ' ->
с очевидными начальными условиями. Границы и центр вихрей определялись из условий иъ (х, г) = 0,
шй(х,г) = 0.
скорости MJ
5м, на линии функции тока ^(х, ^(х)) = Qъ приведены на рис. 3.
При численной реализации метода конечных элементов области бесконечные в горизонтальном
направлении, конечно же, приходится заменять конечными, например при помощи вертикальных границ, имитирующих условие |х| ^ от. Для приведенных результатов расчета выбиралась область —1 < х < 3. На левой границе области х = —1 задавалось условие с = 0, на правой - х = 3. На нижней границе выбирались условия полной проницаемости границ, аналогичные (39), верхней - условие непроницаемости (38).
12 ■ uh u"s\ __. -t / //W X?1^-/i/7/\ f J^' / , % —
, ■ ■ С^ч . ....... / «Л ...... 1, /^Л / "\ ь !/Тх w'\ 6u — !/- x .......
uf 4 ■
—1 w' J? 5и . -1 ■ v^TiY'''/ 3 ---S \ ÖW ■ W
Рис. 2. Придонный вихрь: а - левая граница х = -0,693; z = -0,985; правая граница х = 0,854; z = -1,055; центр вихря х = 0,419; z = -1,115; б - левая граница х = 1,963; z = -0,910; правая граница х = 63,109; z = -0,985; центр вихря х = 2,277; z = -0,956 / Fig. 2. Near-bottom vortex: a -the left boundary is x = -0.693; z = -0.985; the right boundary x = 0.854; z = -1.055; the center of the vortex is x = 0.419; z = -1.115; b - the left boundary is x = 1.963; z = -0.910; the right boundary x = 63.109; z = -0.985; the center of the vortex is x = 2.277; z = -0.956
Зависимости горизонтальной и и вертикальной w скоростей от координаты х на линиях тока, соответствующих показанным на рис. 1, а также средние
Рис. 3. Горизонтальные й£ и вертикальные w£ скорости на линиях функции тока ^£(х,^£) = + Q£ и ^£(х,^) = @й; горизонтальные мй и вертикальные w6 скорости на линиях функции тока ^й(х,^) = и ^£(х,^й + z0) = 0; средние скорости u^v, мйу; величины разрывов скоростей 5u, 5w (40) на линии функции тока ^(х,^(х)) = / Fig. 3. The horizontal u and vertical w velocities on the isolines of the stream function ф'(х,п') = Qb + Q' and ф'(х,г|) = Qb, horizontal ub and vertical wb velocities on the isolines of the stream
function ^b(x,r|) = Qb and ^'(x,r|b + z0) = 0, average velocity u|v, uijv, the magnitude of the velocities discontinuities 5u, 5w (see (40)) on the isoline of the stream function Mx,n(x)) = Qb
Начальное распределение концентрации примеси задавалось в виде пятна, почти полностью локализованного в верхнем слое
c|t=0 = е-^о(х-х°)2, z > - у,
(44)
Ф=о = 0, z<^o-7,
и величины разрывов скоростей
где у > 0, х0 - параметры, характеризующие проникновение примеси в нижний слой, расположение примеси; - параметр сглаживания.
Результаты расчетов положения пятна в различные моменты времени приведены на рис. 4. Для вычислений использовалась триангуляция области, состоящая из 2200 примерно одинаковых треугольников, что соответствует в рассматриваемой области характерному размеру треугольника Ай« 0,04.
ь
и
av
Шаг по времени А С = 0,001 выбирался в соответствии с критерием КФЛ (число Куранта С = (тах(|и|, |ш + з01)А^/А1г и 0,3).
Использование метода конечных элементов с учетом аппроксимации (41) позволяет проводить вычисления с высокой точностью. В частности, закон сохранения массы выполняется с точностью до 1% не только для указанного начального распределения примеси (44), но и практически для любых достаточно гладких распределений.
1.2 1
-10 1 2 3
Рис. 4. Распределение концентрации в различные моменты времени: t = 5 At; t = 85 At-, t = 165 At; t = 255 At; t = 350 At ~ 472,5 s; ß0 = 50; x0 = -0,5; у = 0,01. Внешние границы пятен соответствуют концентрации с = 0,1 / Fig. 4. The concentration distribution at different time: t = 5 At; t = 85 At; t= 165 At; t = 255 At; t = 350 At « « 472.5 s; ß0 = 50; x0 = -0.5; у = 0.01. The spots outer boundary corresponds to the concentration c = 0.1
Заключение
Предложенная модель позволяет достаточно эффективно исследовать поведение пятен примеси в русловых потоках. Может показаться, что модель просто описывает движение слоев жидкости со свободными границами (22), (23), (25). В действительности, помимо того, что движение жидкости в поле тяжести осуществляется за счет уклона потока, выбор вязкости в виде (1), приводящей к сингулярности скорости при z = rjb (x), и использование шероховатости z0(x) позволяют считать границу z = = rb(x) + z0(x) твердой, на которой задается проскальзывание. Во всяком случае, z = rb(x) + z0(x) не является обычной свободной границей жидкости. Наличие явных формул (28), (29) для определения профилей скорости течения и отсутствие ограничений на форму границ слоев дают возможность моделировать практически любую ситуацию. Более того, изучение лишь стационарного течения реально не является препятствием для использования результатов при рассмотрении квазистационарных процессов. Дело в том, что характерные скорости течения
и непосредственно седиментации примеси в поле тяжести, как правило, сильно различаются, и стационарный профиль течения можно с некоторым приближением считать квазистационарным. Следует сказать, что в реальных процессах получить достаточно достоверную информацию, например, на основе измерений можно лишь о рельефе дна водоема (z = ijb(x)) и форме свободной поверхности (z = l(x)). Положение границы раздела между слоями (z = l(x)) вряд ли поддается достоверным измерениям. В некотором смысле это даже преимущество модели, так как анализ задачи показал, что выбор положения границы z = i (x) влияет на наличие либо противотока жидкости, либо придонных вихрей, что, конечно же, уже поддается измерениям или по крайней мере идентификации.
Оправданным является и исследование именно турбулентного течения. Дело в том, что число Рей-нольдса, вычисленное по реальной (не турбулентной) вязкости жидкости, достаточно велико. В приведенных результатах расчетов, если ориентироваться на характерные значения средней скорости течения 0,68 мс-1, глубины 10 м и вязкости 10-6 м2с, число Рейнольдса Re = 6,8 ■ 108 соответствует развитому турбулентному режиму.
Литература
1. Rijn, Leo C. van. Mathematical modeling of morphological processes in the case of suspended sediment transport // Delft Hydra. Commun. 1987. № 382. 260 р.
2. Шлихтинг Г. Теория пограничного слоя. М.: Наука, 1974.
3. ЖуковМ.Ю., ШиряеваЕ.В. Математическое моделирование процесса седиментации примеси в потоке жидкости. Ростов н/Д.: Изд-во ЮФУ, 2016. 208 с.
4. Shiryaeva E.V., Zhukov M.Yu. On the Tangential Stresses at the Boundary Between the Layers for Two-layer Sedimentation Models // Polarforschung. Jahrgang. 2017 (erschienen 2018). Vol. 87, № 2. P. 211-214.
5. Zhukov M.Yu., Shiryaeva E. V. On the completeness problem of the equations for two-layer sedimentation models // Polarforschung. Jahrgang. 2017 (erschienen 2018). Vol. 87, № 2. P. 215-221.
6. Жуков М.Ю., Ширяева Е.В. Влияние седиментации примеси на рельеф дна водоема // Современные проблемы механики сплошной среды: тр. XIX Между-нар. конф.15.10-18.10.2018. Ростов н/Д.; Таганрог: Изд-во ЮФУ, 2018. Т. 2. С. 103-107.
7. Жуков М.Ю., Ширяева Е.В. О касательных напряжениях на границе между слоями для двухслойной модели седиментации примеси в потоке жидкости // Современные проблемы механики сплошной среды: тр. XIX Междунар. конф. 15.10-18.10.2018. Ростов н/Д., Таганрог: Изд-во ЮФУ, 2018. Т. 1. С. 104-108.
8. Бочев М.А., Надолин К.А., Николаев И.А. Моделирование распространения вещества в двумерном стационарном открытом русловом потоке // Мат. моделирование. 1996. Т. 8, № 1. С. 11-24.
9. Монин А.С., Яглом А.М. Статистическая гидромеханика. Механика турбулентности. Ч. I. М.: Наука, 1965. 740 с.
10. Никурадзе И. Закономерности турбулентного движения в гладких трубах // Проблемы турбулентности. М.; Л.: Изд-во ОНТИ НКТП, 1936. С. 75-150.
11. Гришанин К.В. Динамика русловых потоков. Л.: Гидрометеоиздат, 1979. 312 с.
12. Ландау Л.Д., Лифшиц Е.М. Гидродинамика. Теоретическая физика. Т. VI. М.: Наука, 1986. 736 с.
13. CamenenB., LarsonM., BayramA. Equivalent roughness height for plane bed under oscillatory flow // Estuarine, Coastal and Shelf Science. Elsevier, 2009. Vol. 81 (3). 14 p.
14. BelyaevA.V., Vinogradova O.I. Wetting, roughness and flow boundary conditions // J. Phys.: Condens. Matter. 2011. Vol. 23. P. 184104 (15 p.).
15. Lauga E., Brenner M.P., Stone H.A. Microfluidics: The No-Slip Boundary Condition. Chapter 19 // Handbook of Experimental Fluid Dynamics / C. Tropea, A. Yarin, J.F. Foss (Eds.). Springer, 2007. Р. 1219-1240.
16. Жуков М.Ю., Ширяева Е.В. Микрогидродинамика, жидкие плёнки и электрофорез. Ростов н/Д.: Изд-во ЮФУ, 2015. 240 с.
17. Степанянц Ю.А., Фабрикант А.Л. Распространение волн в сдвиговых потоках. М.: Наука, 1996. 240 с.
18. Жуков М.Ю., Ширяева Е.В. Решение задач математической физики при помощи пакета конечных элементов FreeFem++. Ростов н/Д.: Изд-во ЮФУ, 2008. 256 с.
19. Hecht F. New development in FreeFem++ // J. Nu-mer. Math. 2012. Vol. 20, № 3, 4. P. 251-265.
References
1. Rijn, Leo C. van. Mathematical modeling of morphological processes in the case of suspended sediment transport. Delft Hydra. Commun. 1987, No. 382.
2. Shlikhting G. Teoriya pogranichnogo sloya [The theory of the boundary layer]. Moscow: Nauka, 1974.
3. Zhukov M.Yu., Shiryaeva E.V. Matematicheskoe modelirovanie protsessa sedimentatsii primesi v potoke zhidkosti [Mathematical modelling of process of sedimentation of impurities in the fluid stream]. Rostov-on-Don: Izd-vo YuFU, 2016, 208 p.
4. Shiryaeva E.V., Zhukov M.Yu. On the Tangential Stresses at the Boundary Between the Layers for Two-layer Sedimentation Models. Polarforschung. Jahrgang. 2017 (erschienen 2018), vol. 87, No. 2, pp. 211-214.
5. Zhukov M.Yu., Shiryaeva E.V. On the completeness problem of the equations for two-layer sedimentation models. Polarforschung. Jahrgang. 2017 (erschienen 2018), vol. 87, No. 2, pp. 215-221.
Поступила в редакцию /Received_
6. Zhukov M.Yu., Shiryaeva E.V. [The effect of sedimentation of impurities on the bottom relief of the reservoir]. Sovremennye problemy mekhaniki sploshnoi sredy [Modern problems of continuum mechanics]. Proceedings of the 19th International Conference 15.10-18.10.2018. Rostov-on-Don; Taganrog: Izd-vo YuFU, 2018, vol. 2, pp. 103-107.
7. Zhukov M.Yu., Shiryaeva E.V. [On shear stresses at the interface between layers for a two-layer model of impurity sedimentation in a fluid flow]. Sovremennye problemy mekhaniki sploshnoi sredy [Modern problems of continuum mechanics]. Proceedings of the 19th International Conference 15.10-18.10.2018. Rostov-on-Don; Taganrog: Izd-vo YuFU, 2018, vol. 1, pp. 104-108.
8. Bochev M.A., Nadolin K.A., Nikolaev I.A. Mod-elirovanie rasprostraneniya veshchestva v dvumernom statsionarnom otkrytom ruslovom potoke [Modeling of substance propagation in a two-dimensional stationary open channel flow]. Mat. modelirovanie. 1996, vol. 8, No. 1, pp. 11-24.
9. Monin A.S., Yaglom A.M. Statisticheskayagidromek-hanika. Mekhanika turbulentnosti [Statistical hydromechanics. Mechanics of turbulence]. Ch. I. Moscow: Nauka, 1965.
10. Nikuradze I. [Regularities of turbulent motion in smooth pipes]. Problemy turbulentnosti [Turbulence problems]. Moscow; Leningrad: Izd-vo ONTI NKTP, 1936, pp. 75150.
11. Grishanin K.V. Dinamika ruslovykh potokov [Dynamics of channel flows]. Leningrad: Gidrometeoizdat, 1979.
12. Landau L.D., Lifshits E.M. Gidrodinamika. Teoret-icheskaya fizika [Hydrodynamics. Theoretical physics]. Vol. VI. Moscow: Nauka, 1986.
13. Camenen B., Larson M., Bayram A. Equivalent roughness height for plane bed under oscillatory flow. Estuarine, Coastal and Shelf Science. Elsevier, 2009, vol. 81 (3), 14 p.
14. Belyaev A.V., Vinogradova O.I. Wetting, roughness and flow boundary conditions. J. Phys.: Condens. Matter. 2011, vol. 23, p. 184104 (15 p.).
15. Lauga E., Brenner M.P., Stone H.A. Microfluidics: The No-Slip Boundary Condition. Chapter 19. Handbook of Experimental Fluid Dynamics. C. Tropea, A. Yarin, J. F. Foss (Eds.). Springer, 2007.
16. Zhukov M.Yu., Shiryaeva E.V. Mikrogidro-dinamika, zhidkie plenki i elektroforez [Microhydrody-namics, liquid films and electrophoresis]. Rostov-on-Don: Izd-vo YuFU, 2015, 240 p.
17. Stepanyants Yu.A., Fabrikant A.L. Raspros-tranenie voln v sdvigovykh potokakh [Wave propagation in shear flows]. Moscow: Nauka, 1996, 240 p.
18. Zhukov M.Yu., Shiryaeva E.V. Reshenie zadach matematicheskoi fiziki pri pomoshchi paketa konechnykh elementov FreeFem++ [Solving problems of mathematical physics using the finite element package FreeFem++]. Rostov-on-Don: Izd-vo YuFU, 2008, 256 p.
19. Hecht F. New development in FreeFem++. J. Nu-mer. Math. 2012, vol. 20, No. 3, 4, pp. 251-265.
_4 июня 2019 г. / June 4, 2019