Self-similar solutions of the bed deformation

全文:

Традиционно для изучения деформаций дна, возникающих под воздействием протекающего по нему гидродинамического потока, используют нестационарные русловые модели [1], требующие значительных вычислительных затрат [2]. Поскольку результаты решения таких задач используются в большом количестве прикладных проектов (прокладка трубопроводов по дну рек [3, 4], деформирование донной поверхности в нижних бьефах гидроузлов [5–7], воздействия движителей судов на дно [8] и т.д.), то возможность сокращения времени расчета имеет большое значение для практики.

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

1. Автомодельный размыв дна

Экспериментальные данные донной поверхности, размываемой под действием гидродинамического потока [3–8], показывают, что данные процессы имеют автомодельный характер [9–11]. На рис. 1 представлены серии экспериментальных кривых, из которых видно, что донную поверхность во все моменты времени можно представить в следующем виде:

$\sigma \left(\frac{x}{{\zeta }_0\left(t\right)}\right)=\frac{\zeta (t,x)}{{\zeta }_0\left(t\right)}=\left\{\begin{array}{c}-\sin 2\pi q\frac{x}{{\zeta }_0\left(t\right)},q\frac{x}{{\zeta }_0\left(t\right)}\in \left(\mathrm{0,}1\right)\\ \mathrm{0,}q\frac{x}{{\zeta }_0\left(t\right)}\notin \left(\mathrm{0,}1\right)\end{array}\right.$ (1)

$q={\zeta }_0\text{/}\lambda .$ (2)

Физически величина q равна крутизне донной волны и по данным экспериментов различных русловых процессов является постоянной величиной. Например, по данным экспериментальных кривых для эволюции донной поверхности под трубой на дне реки q = ζ₀ / λ = 0.08 (рис. 1а [11]). Другой пример эволюции донной поверхности под воздействием плоской придонной струи (течение из-под щита), по данным работы [6], изображен на рис. 1б. Из него видно, что для шести различных моментов времени сохраняется синусоидальная форма донной поверхности, меняются длина и амплитуда донной волны, а их отношение остается постоянным.

Рис. 1. Эволюция донной поверхности: а – под трубой на дне реки [11], б – под воздействием плоской придонной струи по данным работы [6], 1 – синусоидальная аппроксимация.

Автомодельные синусоидальные формы (1) изображены на рис. 1 в виде сплошных линий. Изменение амплитуды волны ζ₀(t) и ее длины λ(t) во времени аппроксимируются степенными функциями

$\lambda =\alpha t^{\beta },{\zeta }_0=q\alpha t^{\beta }.$ (3)

Например, для автомодельных донных форм на рис. 1б α = 70, β = 1 / 8.

2. Математическая постановка задачи

При описании эволюции уровня донной поверхности ζ используется закон сохранения донного материала (уравнение Экснера)

$\left(1-\epsilon \right)\frac{\partial \zeta }{\partial t}+\frac{\partial G}{\partial x}=\mathrm{0,}$ (4)

где G – объемный расход влекомых наносов для частиц грунта, движущихся в активном придонном слое, ε – пористость донного материала.

Для замыкания уравнения (4) используется теоретическая зависимость

$\begin{array}{c}G=G_0\frac{{\left(T-T_c\right)}^{3/2}\eta \left(T-T_c\right)}{1+\Gamma },\Gamma =\frac{1}{tg\varphi }\frac{\partial \zeta }{\partial x},\\ \\ T_c=T_{0c}\left(1+\Gamma \right),\\ G_0=\frac{g_0}{g\left({\rho }_s-{\rho }_w\right)\sqrt{{\rho }_w}},g_0=\frac{4}{3\kappa tg\phi },\\ \\ T_{0c}=\frac{9}{8}\frac{{\kappa }^{2}tg\varphi \left({\rho }_s-{\rho }_w\right)gd}{c_x},\end{array}$ (5)

где T_c – критическое придонное касательное напряжение, T_c0 – критическое напряжение на ровном дне, η(z) – функция Хевисайда.

Теоретическая зависимость (5) подтверждена всеми наиболее известными и надежными экспериментальными данными. Уравнения (5) и (4) совместно с уравнениями гидродинамики представляют собой замкнутую, но весьма сложную систему уравнений, позволяющую рассчитывать эволюцию донной поверхности под воздействием различных гидродинамических факторов: течение под трубой, под воздействием придонной струи и другие. Используя автомодельную геометрию донной волны, можно описать ее эволюцию, сократив расчетное время почти на два порядка.

3. Уравнение для распределения напряжения на донной поверхности при отсутствии транзитных наносов

С помощью автомодельного решения (1) можно найти распределение напряжения на донной поверхности. Для этого подставим его в уравнение Экснера. Для производной по времени функции ζ(t, x) имеем

$\begin{array}{l}\frac{\partial \zeta }{\partial t}=q\left(-\dot{\lambda }\sin \left(\frac{2\pi }{\lambda \left(t\right)}x\right)+\lambda \cos \left(\frac{2\pi }{\lambda \left(t\right)}x\right)\frac{2\pi x}{{\lambda }^2\left(t\right)}\dot{\lambda }\right)=\\ =q\dot{\lambda }\left(-\sin \left(\frac{2\pi }{\lambda \left(t\right)}x\right)+\frac{2\pi x}{\lambda \left(t\right)}\cos \left(\frac{2\pi }{\lambda \left(t\right)}x\right)\right),\dot{\lambda }=\frac{d\lambda }{dt}.\end{array}$

Проинтегрируем уравнение Экснера (4) по x:

$\begin{array}{c}(1-\epsilon )\underset{0}{\overset{x}{\int }}\frac{\partial \zeta }{\partial t}dx+G(t,x)=G\left(0\right),\\ \underset{0}{\overset{x}{\int }}\frac{\partial \zeta }{\partial t}dx=q\dot{\lambda }\left[x\sin \left(\frac{2\pi x}{\lambda }\right)-\frac{\lambda \left(1-\cos \left(\frac{2\pi x}{\lambda }\right)\right)}{\pi }\right]=\\ =-\frac{q}{2\pi }\dot{\lambda }\lambda \Phi \left(\xi \right)\Rightarrow \\ \Rightarrow G(t,x)-G\left(0\right)=(1-\epsilon )\frac{q}{2\pi }\dot{\lambda }\left(t\right)\lambda \left(t\right)\Phi \left(\xi \right),\\ \Phi \left(\xi \right)=2-2\cos \xi -\xi \sin \xi ,\xi =\frac{2\pi x}{\lambda }.\end{array}$

Отсюда, предполагая, что G(0) – (транзитный расход наносов) равен нулю, с помощью (5) получим уравнение для распределения напряжения на донной поверхности T:

$G_0\frac{{\left(T-T_c\right)}^{3/2}\eta \left(T-T_c\right)}{1+\Gamma }=(1-\epsilon )\frac{q}{2\pi }\dot{\lambda }\left(t\right)\lambda \left(t\right)\Phi \left(\xi \right),$

разрешая которое относительно T, найдем

$\begin{array}{c}\frac{T}{T_{0c}}=1-\frac{2\pi q}{tg\phi }\cos \xi +\\ +{\left((1+\Gamma (\xi \left)\right)\Phi \left(\xi \right)\right)}^{2/3} \Lambda {\left(t\right)}^{2/3},\\ \Lambda \left(t\right)=\frac{\left(1-\epsilon \right)q}{2\pi }\frac{\dot{\lambda }\left(t\right)\lambda \left(t\right)}{G_0{T_{0c}}^{3/2}},\\ \Gamma =\frac{1}{tg\varphi }\frac{\partial \zeta }{\partial x}=-\frac{2\pi q}{tg\phi }\cos \xi ,\\ \xi =\left(\frac{2\pi }{\lambda \left(t\right)}x\right).\end{array}$ (6)

Графики функций Φ(ζ) и Г(ζ), входящих в (6), изображены на рис. 2а, 2б соответственно.

Рис. 2. Графики функций Ф(ξ) (а) и Г(ξ) (б).

С помощью (6) можно получить связь среднего значения напряжения на донной волне с функцией Λ(t):

$\begin{array}{l}\frac{\overline{T}\left(\lambda \right)}{T_{0c}}=\frac{1}{2\pi }\underset{0}{\overset{2\pi }{\int }}\frac{T}{T_{0c}}d\xi =1+-{\Lambda }^{2/3},\\ -=\frac{1}{2\pi }\underset{0}{\overset{2\pi }{\int }}{\left(1-\frac{2\pi q}{tg\phi }\cos \xi \right)}^{2/3}\Phi {\left(\xi \right)}^{2/3}d\xi .\end{array}$

Коэффициент С при изменении множителя $\frac{2\pi q}{tg\phi }$ на интервале (0.4, 0.8) меняется незначительно в пределах (1.97, 2.01). Поэтому можно принять

$\frac{\overline{T}\left(\lambda \right)}{T_{0c}}=\frac{1}{2\pi }\underset{0}{\overset{2\pi }{\int }}\frac{T}{T_{0c}}d\xi =1+2{\Lambda }^{2/3}.$ (7)

Разрешая с помощью (6) это уравнение относительно dλ / dt, получим дифференциальное соотношение:

$\frac{d\lambda }{dt}=\frac{K}{\lambda }{\left(\frac{\overline{T}\left(\lambda \right)}{T_{0c}}-1\right)}^{3/2},K=\frac{\pi G_0{T_{0c}}^{3/2}}{\sqrt{2}(1-\epsilon )q}.$ (8)

4. Схема расчета эволюции автомодельной донной волны

Предлагается следующий алгоритм для расчета эволюции донной волны в автомодельной форме (3) с помощью уравнения (8). Пусть имеется набор донных волн в виде (1) с длинами волн λ₁, λ₂, ..., λ_k. Для каждой донной волны λ_i из решения уравнений гидродинамики численно находится распределение напряжений на равномерной сетке донной волны: T₁, T₂, ..., T_N. По этим значениям по формулам

$\frac{\overline{T}\left(\lambda \right)}{T_{0c}}=\frac{1}{N}\sum _{i=1}^N\frac{T_i}{T_{0c}},T_{0c}=\frac{9{\kappa }^{2}tg\varphi }{8c_x}({\rho }_s-{\rho }_w)gd$

вычисляются значения

$\overline{T}\left({\lambda }_1\right)/T_{0c},\overline{T}\left({\lambda }_2\right)/T_{0c},\mathrm{...,}\overline{T}\left({\lambda }_k\right)/T_{0c}$

и из дифференциального соотношения (8) находится таблица значений производных

${\dot{\lambda }}_1\left({\lambda }_1\right),{\dot{\lambda }}_2\left({\lambda }_2\right),\mathrm{...,}{\dot{\lambda }}_k\left({\lambda }_k\right),{\dot{\lambda }}_i=\frac{d{\lambda }_i}{dt}.$

Аппроксимация табличных значений в виде

$\frac{d\lambda }{dt}=A{\lambda }^{-n}$ (9)

находится методом наименьших квадратов. Из условия минимума функции

$F=\sum _{i=1}^k(\ln {\dot{\lambda }}_i-\ln A+n\ln {\lambda }_i{)}^2$

получаем систему двух линейных относительно lnA и n уравнений

$\begin{array}{c}\sum _{i=1}^k(\ln {\dot{\lambda }}_i-\ln A+n\ln {\lambda }_i)=\mathrm{0,}\\ \sum _{i=1}^k(\ln {\dot{\lambda }}_i-\ln A+n\ln {\lambda }_i)\ln {\lambda }_i=0.\end{array}$  (10)

В свою очередь, уравнение (9) имеет точное решение

$\lambda \left(t\right)={\left(\left(n+1\right)A\right)}^{1/\left(n+1\right)}{t}^{1/\left(n+1\right)}.$  (11)

После подстановки в него найденных из системы уравнений (10) значений A и n получаем функцию λ(t). Такой подход позволяет существенно сократить расчетное время.

5. Расчетный пример

Для численного расчета напряжений на синусоидальном дне использовался открытый программный комплекс OpenFOAM. Рассматривались обтекания цилиндра (диаметр D = 0.1 м) с различными формами донной поверхности потоком воды со скоростью 0.4 м/c. Определялись напряжения на донной поверхности для синусоидальных форм с длинами волн λ = 0.125, 0.2, 0.3, 0.4, 0.6 и крутизной q = a / λ = 0.08. На рис. 3 представлены постановка задачи (а) и различные формы донных синусоидальных волн (б); результаты расчета обтекания цилиндра над донными волнами: длина волны 3D (в) и длина волны 6D (г).

Рис. 3. Постановка задачи обтекания цилиндра: а – схема расчетной области. Г_in – входная граница, Г_out – выходная граница, Г_bed – донная граница, Г_top – верхняя граница; б – формы донной поверхности; результаты расчета обтекания цилиндра над донными волнами: в – длина волны 3D; г – длина волны 6D.

По результатам численного моделирования для пяти значений длин волн λ с крутизной волн q = 0.08 получены средние напряжения $\overline{T}=\mathrm{0.99,}\mathrm{0.84,}\mathrm{0.56,}\mathrm{0.46,}0.42$. Данные для песка на дне в системе СИ: диаметр песчинок d = 0.00035, пористость ε = 0.35, параметр Кармана ê = 0.25, коэффициент сопротивления песчаных частиц c_x = 0.55, предельный угол трения tg φ = 0.5, плотность песка ρ_s = 2650, плотность воды ρ_w = 1000. Напряжение на бесконечности ~0.36. Транзитные наносы отсутствуют.

По представленным данным определяем

$\begin{array}{c}{G}_0=\frac{4}{3\kappa tg\phi g\left({\rho }_s-{\rho }_w\right)\sqrt{{\rho }_w}}=\mathrm{0.0000209,}\\ T_{0c}=\frac{9{\kappa }^{2}tg\varphi }{8c_x}({\rho }_s-{\rho }_w)gd=\mathrm{0.362,}\\ K=\frac{\pi G_0{T_{0c}}^{3/2}}{\sqrt{2}(1-\epsilon )q}=0.000194\end{array}$

и с помощью (8) находим значения скорости изменения длины волны dλ / dt = 0.00355, 0.00147, 0.000262, 0.0000686, 0.0000209.

Результаты расчетов представлены в табл. 1.

Таблица 1

Методом наименьших квадратов из системы уравнений (10) находим коэффициенты интерполяционной функции $\dot{\lambda }=A{\lambda }^{-n}$ lnA = −12.47, A = 0.00000385, n = 3.42. Сравнение интерполяционной кривой $\dot{\lambda }=A{\lambda }^{-n}$ с представленными в табл. 1 красными точками приведено на рис. 4а. Решение определяется по формуле (11). Подставляя найденные значения A, n, получим

$\lambda \left(t\right)\left(м\right)\text{\hspace{0.05em} =\hspace{0.05em} }0.0834{\left(t\left(c\right)\right)}^{0.226}\text{ = \hspace{0.17em}}0.211{\left(t\left(мин\right)\right)}^{0.226}$ .

График зависимости

$\zeta \left(t\right)/D=0.0169t^{0.226}/D$

изображен на рис. 4б синей линией в сравнении с экспериментальными данными [11] и результатами расчета непрерывного размыва донной поверхности и деформации сетки [12].

Рис. 4. Сравнение результатов: а – расчетные точки со степенной аппроксимацией; б – расчет [12] (1), экспериментальные данные [11] (красные точки), предложенная модель (2).

Из графика видно удовлетворительное совпадение численного расчета и эксперимента с теоретической моделью. Погрешность вычислений ζ(t) не превышает нескольких процентов. При этом численное моделирование гидродинамики с учетом размыва требует до 8 ч. Определение напряжений для 5 рассмотренных случаев и решение задачи по предложенному алгоритму составило порядка 15 мин, т.е. расчетное время сократилось в 30 раз. Стоит отметить, что предложенная модель зависит от численных результатов напряжений на донной поверхности, которые напрямую определяются моделью турбулентности и их коэффициентами, реализацией пристеночных функций, расчетной сеткой, что может привести к увеличению погрешности расчета. Можно также отметить, что точность расчета практически не изменится, если выбрать из таблицы значения только для четырех длин волн $\lambda =:1.25D,3D,4D$. При этом время сократится на 20%.

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

Источник финансирования

Исследования выполнены за счет гранта Российского научного фонда № 23-71-10091.

作者简介

А. Petrov

编辑信件的主要联系方式.
Email: petrovipmech@gmail.com
俄罗斯联邦, Moscow

I. Potapov

Email: petrovipmech@gmail.com

Хабаровский федеральный исследовательский центр

A. Epikhin

Email: petrovipmech@gmail.com
俄罗斯联邦, Moscow

参考

补充文件

附件文件

动作

1. JATS XML

下载

2. Fig. 1. Evolution of the bottom surface: a – under a pipe on the river bottom [11], b – under the influence of a flat bottom jet according to data from [6], 1 – sinusoidal approximation.

下载 (191KB)

索引源数据

3. Fig. 2. Graphs of functions Ф(ξ) (a) and Г(ξ) (b).

下载 (99KB)

索引源数据

4. Fig. 3. Statement of the problem of flow past a cylinder: a – scheme of the computational domain. Гin – input boundary, Гout – output boundary, Гbed – bottom boundary, Гtop – upper boundary; b – shapes of the bottom surface; results of calculation of flow past a cylinder over bottom waves: c – 3D wavelength; g – 6D wavelength.

下载 (388KB)

索引源数据

5. Fig. 4. Comparison of results: a – calculated points with power approximation; b – calculation [12] (1), experimental data [11] (red dots), proposed model (2).

下载 (201KB)

索引源数据

Presented by Academician of the RAS R.I. Nigmatullin