CC BY
22
УДК 551.345 + 539.3
Численное решение двумерной задачи фильтрации в верхних слоях почвогрунтов с учетом суффозионных процессов *
А.Н. Сибин, Н.Н. Сибин
Алтайский государственный университет (Барнаул, Россия)
Numerical Solution of Two-Dimensional Filtration Problem in the Upper Layers of Soil Considering Suffosion Processes
A.N. Sibin, N.N. Sibin
Altai State University (Barnaul, Russia)
Рассматривается математическая модель изотермической внутренней эрозии без учета деформации пористой среды. При достижении определенной величины скорости фильтрации происходит вынос частиц грунта из области течения. В качестве математической модели используются уравнения сохранения массы для воды, подвижных твердых частиц и неподвижного пористого скелета, а также закон Дарси для воды и подвижных твердых частиц и соотношение для интенсивности суффозионного потока. В пункте 1 дается постановка задачи и проводится преобразование системы уравнений. В результате преобразований для насыщенности водной фазы возникает вырождающееся на решении параболическое уравнение для давления эллиптическое уравнение и уравнение первого порядка для пористости грунта. Имеется аналогия с классической моделью Маскета - Леверетта. В пункте 2 предложен алгоритм численного решения двумерной начально-краевой задачи внутренней эрозии грунта. В пункте 3 представлены результаты численного решения задачи. Найдена область, наиболее подверженная внутренней суффозии.
Ключевые слова: многофазная фильтрация, пористая среда, суффозия, фазовый переход, насыщенность.
БО! 10.14258^уа8и(2017)4-27
In this paper mathematical model of isothermal internal erosion without deformation of porous medium is studied. Removal of soil particles from flow occurs at a certain magnitude of velocity of filtration. Equations of mass conservation of water, moving solids and stationary porous skeleton along with Darcy's law for water and moving solid particles and equation for the intensity of suffusion aquifer were used as the mathematical model of the problem. Formulation of the problem and development of the system of equations are described in section 1. The results are parabolic equation with extinction of solution, elliptical equation for pressure and equation of the first order for the porosity of the soil. There is analogy with the classical Masket-Leverett model. Algorithm of numerical solution of two-dimensional initial boundary value problem of internal soil erosion proposed in section 2. The results of numerical investigation of the problem presented in section 3. We found the region of most susceptible to internal suffosion.
Key words: multiphase flow, porous medium, suffusion, phase transition, saturation.
1. Постановка задачи. Проблемы, связанные с эрозией почвогрунтов, широко исследуются на протяжении последнего столетия. Данный процесс имеет большое значение при решении прикладных задач в сельском хозяйстве: ирригация и дренаж сельскохозяйственных полей [1] и процесс внутренней эрозии, сопутствующий канальному орошению почвогрунтов [2]. Процесс
* Работа выполнена при финансовой поддержке грантов РФФИ №16-08-00291 и №17-41-220314.
эрозии необходимо учитывать в исследованиях, связанных с прогнозом распространения загрязнений, фильтрацией вблизи водохранилищ и других гидротехнических сооружений [3-5]. Более того, аналогичные проблемы, связанные с процессом эрозии грунта, возникают и в других областях, включая добычу нефти и газа [6].
Эрозия почв отрицательно влияет практически на все отрасли сельскохозяйственного производства. Наряду со снижением плодородия почв,
Численное решение двумерной задачи фильтрации.
уменьшением урожая культур и другими вредными воздействиями на сельскохозяйственное производство эрозия значительно ухудшает условия функционирования сельскохозяйственных машин и агрегатов при выполнении ими технологических работ. Это проявляется в снижении производительности работы сельскохозяйственной техники на склонах, расчлененных промоинами и оврагами, в ухудшении качества полевых работ, увеличении износа машин [1,2].
В данной работе рассматриваются процессы фильтрации воды в верхних слоях почвогрунтов и внутренней суффозии. Почвогрунт моделируется как трехфазная сплошная пористая среда. Поры полностью заполнены смесью воды (1 = 1) и подвижных твердых частиц (1 = 2). Доля пор в грунте (1 = 3) определяется пористостью ф = (VI + VI)¡V, где V = VI + V2 + У3 - общий объем грунта, Vl, V2, Уз - соответственно объемы воды, подвижных твердых частиц и скелета грунта.
Следуя работе [7], в двумерном случае обезраз-меренная система уравнений, описывающая процесс внутренней эрозии, имеет вид
дз
ф— = У-(аУз + Ъу + F),
У-(К Ур + /)=0,
дф = т дь т,
(1)
(2)
(3)
где I - интенсивность фазового перехода (суффо-зионный поток); .з = V1¡(V1 + V2), .з2 = V2¡(V1 + ^2) - концентрации воды (насыщенность) и подвижных твердых частиц в порах; Рх,Р2,Рз - истинные плотности воды, подвижных твердых частиц грунта и скелета грунта; V = (дх, ду) - оператор градиента; V = ("х,"у) - суммарная скорость фильтрации, р - приведенное давление. В рассматриваемом случае р° = р!°, так как подвижные частицы захватываются суффозионным потоком из грунта. Коэффициенты уравнения имеют вид:
^ к02 дР2 р2
""у2 = -Кот (дх - рзз ду);
к01 Ы р - 1)
^ = Кодх-у—-;
Тк01 + уз2
Кз
В р доЪ 3
—2— V , ^ = Кзду
хвси1
р2Т1
к01 к02 ( % - 1) Тк01 + к02
¡х = -Ко(уо1 + -„— кз2)дх, К = Кок, Р1Т 2
к = у01 + — У02, ¡у = -Ко(уо1 + ко2)ду.
Т2 Р1Т 2
А(1 - ф)(1 - .з)ф(\"1\-^ |), \v1\y\vk |;
0, Ы < \vk\.
А = Ах6с, vk = vk —, \"1 \ = ^ +
0 2
- ^ V = х
размер-
Здесь Ьвс,Хвс,рвс ¿2 , ивс
ные постоянные; АА - размерный эмпирический коэффициент (отвечает за суффозионную устойчивость почвогрунта), а ц = ^.
2. Алгоритм численного решения задачи. Введем сетку с распределенными узлами хг = ihl, yj = jh2, Ьп = пт; i = 0,..., N1, j = 0,..., N2 п = 0,...,Т, Ь, h2
- шаг по пространственным координатам, т
- шаг по времени, итерационный параметр
= ^5аЛа, 5а = -щ sin2(Па), Ла =
-42 cos2(Па), а = 1, 2; 11 и 12 длина и ширина области соответственно, N1, ^-число разбиений 11 и ¡2.
Уравнения (1) и (2) решаем методом переменных направлений [8]. Разностная схема имеет следующий вид:
^1Рп+1 - Л1Рп+1 = ^1Рп - Л2Рп + Ф; (4)
^2Рп+1 - Л2Рп+1 = ^2Рп+1 - Л1Рп+1 + ф; (5)
, Ч ^ к01у02 др^ к02
Тк01 + к02 дв Тк01 + к02
х вс . х вс . Л1Р
"х = — = + "х2, Vу = Vу — = "у1 + "у2 ;
ьвс ьвс
"х1 = -Коко1 (дХ - дх);
к02 др2 р02 "х2 = -Ко-(^---о дх );
дх р1
"у1 = Коко1(дХ - ду);
п+ 2 п+ 2 п+ 2 п+ 2
= кп р+л - ры _ К п Ры - р-М ;
, Н\ 2 л h2 :
1
пп
пп
^^ _^ ^ _р
2 гл+2 h2 г'3-1 h2
Ф = ¡хфг,.
кп _ Хп оп _ оп
2h1
+ ¡хвг,л
2h1
+
т
вс
dF рп,. . = хф^3 дф (вп■ фп •) дР Рп =W1 y ( п фп ). y,i,j дф ,j i,jh
f п dfx fx,i,j = дф (вп■ фп ■) f п = д fy (вп фп ). Jy,i,j дф i,jlVi,jh
fn dfx Jxsi'j = дв (вп■ фп ■) fn = dfy(вп ф" ). Jysi,J дв ,j j
(6) v ' Gn Gxi,j dF —(вп. ф" дЬ j )+ Vx(tn) дГз (в" У.
Gn Gyi,j = dF (вп. ф" дЬ j )+ Vy (t") g~s (в" У.
K" 2 • = 2 j 2K (ф"±и )K (фпj )
K (фп±и ,в n±1,j) + K (фп ).
K" = Knj± 2 = 2K (ф"± , вп^±1)к(фп , )
K (ф"л±1,в "j±1) + K (фпj ,впj).
+fyфi,:
,п - фп
i,j + 1
i,j-1
2h
+ fy
пп si,j + 1 si,j-1
2
ysi,j
2h
2
При аппроксимации уравнения (1) используется направленная разность для конвективного слагаемого.
п+2
в- ■ — в- ■ ФЬ i,j T ij = Г 1вп+ 2 + Г2вп + п - в.
фп+1 фп ij T
вп+1 _ „п+ 2
in Si,j Si,j T
+п - в.
= Г2 вп+1 +Г1 вп+2 + фп+1 _ фп .
(7)
Г1в
п+2 п+2
о 2 - Q 2
п si+1,j si,j ап a - | i --гт;- — a - i -
i+2 j h2 i-1 j
п+2 п+ 2
Q 2 - Q 2
si,j si-1,j
hi
+
+
¡. п+ 2 r> |/On I п+ 2 i Л п+ 2 i1sei+1,j — 2\Gxi,j\si,j + >2x si-1,j
2h1 '
п п п п
> 1 x = , \Gxi,j \ + Gxi,j S2x = \Gxi,j \ - G
xi,j
Г2 в
пп ап bi,j + 1 bi,j ап a- - i 2-гт;- — a- - 2
1J + 2 h2 i,j-2
qio _ qio
bi,j 1
h2
+
+
^ в".+1 - 2(\G"j Mjj + by j 2h2 '
п п п п
s1y = \Gyi,j \ + Gyi, j , S2y = \Gyi,j \ Gyi,
п = fx
Фп+Uj - фп
i- 1 ,j
фnj+l - фпj
+ ■-
2h1 + y,ij 2h2
i,j-1
Уравнение 3 аппроксимируется (неявной) схемой Рунге - Кутта второго порядка точности. Причем найденное на первом этапе значение фП+1
ФП+1 = фЬ + Tin
затем уточняется следующим образом:
I №j ) +1 (ф>п< )
фп+1 = фЬ + T-
2
(8)
(9)
Здесь i = 1,..,N - 1j = 1,..,N2 -1, t =1,..,T -1,
2Ко(ф±и )а(вП±и )Ko№j )а(вП<)
± 2 ,j Ко(фП±и )а(вП±и) + Ко(фПj )а(впj)
о(фПл )а(вЬ>
А(1 - ф^ )(1 - П )фп (К, \-\vu |),
К,\ \;
о, К) \ < К \.
= , ) )2 + , ) ?.
Алгоритм численного решения начально-краевой задачи следующий: используя начальное значение пористости ф0 ^ и концентрации в0) находим начальное распределение приведенного давления р0 ^ (г = 0,...,М\;j = 0,...,М2), решив систему разностных уравнений (4) и (5) методом переменных направлений. Используя найденное давление, определяем скорости фильтрации воды и подвижных частиц грунта ) и ). Из равенства (8) находим пористость грунта ф1) на следующем шаге по времени. Находим концентрацию воды в1), решив систему разностных уравнений (6) и (7) методом переменных направлений, используя найденные давления и пористость.
ai, j± 2
2Ко(фП±)а(вП±)Ко(фП^ )а(впj) 2 Ко(фП±)а(вП±) + Ко(фПj)а(впj)'
Изменение пористости грунта
п
Численное решение двумерной задачи фильтрации...
Рассчитываем давление на следующем шаге по времени (п = 1). Используя найденные значения искомых функций ф1 ^, я1 ^, р0 ^, р1 ^, корректируем значения пористости на первом шаге по времени, используя формулу (9). Повторяя данный алгоритм для следующих шагов по времени, найдем значения искомых функций на всем временном интервале.
3. Результаты численных расчетов. На поверхносте грунта поддерживается слой воды глубиной Лводы, достаточной для возникновения суффозионного процесса. Предпологается, что на нижней границе рассматриваемой области суфо-зионный процесс не происходит, так как скорость воды меньше критической скорости, необходимой для возникновения суффозионного процесса.
Начальные и граничные условия задачи имеют вид
дя дя
s(x,y, 0) = 1, — (0,у,*) = 0 — (1ьу,*)=0;
дХ(х, 0,^=0 я(х,М) = 1;
др(0,у,*) = 0 др(11,у,*) = 0;
дХХ(х,0,^ = 0 р(х,г2,г) = Л-водыРоя;
ф(ху ») = {аз,хх> !:воды:
В численных расчетах использовался следующий набор модельных параметров [9]:
Р0 = 1000 кг/м3; р0 = 1440 кг/м3;
В = 0:00001 м2; я = 9:8 м/с2, ук = 0;
= 0:001787 кг/мс; ц2 = 0:003574 кг/мс
На рисунке представлен график изменения пористости грунта при х = 0:5. С течением времени выделяется область, где почва размывается полностью (ф = 1) (см. рис.). Численное решение одномерной задачи фильтрации грунтовых вод с учетом суффозионных процессов получено в работе [10].
Заключение. В работе дан краткий литературный обзор моделей внутренней эрозии почвы. Предложен алгоритм численного решения двумерной задачи фильтрации воды в верхних слоях почвогрунта с учетом внутренней суффозии и проведены тестовые численные расчеты. Найдена область, наиболее подверженная внутренней суффозии.
Авторы статьи признательны А.А. Папину и С.С. Кузикову за обсуждение задачи и конструктивные замечания.
Библиографический список
1. Vieira D.A.N., Dabney S.M. Modeling edge effects of tillage erosion // Soil Tillage Research. — 2011. 111(2).
2. Wilson G. Understanding soil-pipe flow and its role in ephemeral gully erosion // Hydrol. Process. — 2011. — Vol. 25.
3. Einstein H. A. Der Geschiebetrieb als wahrscheinlichkeits Problem. Mitt. d. Versuchsanstaltf Wasserbau, Eidg. T. H. — Zurich, 1937.
4. Bonelli S. Erosion of Geomaterials. UK, 2012.
5. Vardoulakis I., Stavropoulou M., Papanastasiou R. Hydro-Mechanical Aspects of the Sand Production Problem // Transport in Porous Media. — 1996. — 22.
6. Wang J., Walters D. A., Settari A., Wan R. G. Simulation of cold heavy oil production using
an integrated modular approach with emphasis on foamy oil flow and sand production effects 1st Heavy Oil Conference. - 2006.
7. Papin A. A., Sibin A. N. Model isothermal internal erosion of soil //J. Phys.: Conf. Ser. — 2016. — V. 722(1).
8. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. — М., 1978.
9. Chetti A., Benamar A., Hazzab A. Modeling of Particle Migration in Porous Media:Application to Soil Suffusion // Transport in Porous Media. — 2016. — V. 113(3).
10. Сибин А.Н., Сибин Н.Н. Расчет физических характеристик почвогрунтов в процессе внутренней эрозии // МАК: Математики — Алтайскому краю : сборник трудов Всероссийской конференции по математике. — Барнаул, 2017.
