CC BY
46
Мурзакматов М.У.
Кафедра прикладной математики Иссык-Кульского государственного университета,
кандидат технических наук, доцент Исабеков К.А.
Кафедра прикладной математики Иссык-Кульского государственного университета
ПРИБЛИЖЕННОЕ РЕШЕНИЕ КРАЕВОЙ ЗАДАЧИ ПРОСТРАНСТВЕННОЙ СТАЦИОНАРНОЙ ФИЛЬТРАЦИИ ПОДЗЕМНЫХ ВОД
Рассматривается вывод уравнения стационарного пространственного течения подземных вод, основанный на гидродинамическом уравнении Эйлера и законе Дарси, и разрабатывается алгоритм его приближенного решения методом конечных элементов. Алгоритм и программа отлажены на решении тестовых задач.
При выводе уравнений движения подземных вод исходят из следующих упрощающих допущений: жидкость, полностью заполняющая поры грунта в области течения, однородна и несжимаема; под водоносной средой понимается грунт, отдельные частицы которого неподвижны и устойчивы в своем взаимном расположении, при этом грунт считается несжимаемым; капиллярные силы вдоль свободной поверхности не учитываются; физические параметры водоносного пласта (пористость, водо-проводимость, коэффициент фильтрации и др.) считаются не изменяющимися с течением времени; силы сопротивления зависят только от трения частиц жидкости о частицы грунта; и, наконец, допускается, что скорости фильтрации невелики, силы инерции пренебрежимо малы по сравнению с силой тяжести и силами трения.
При этих предложениях гидродинамическое уравнение Эйлера имеет вид:
(1)
где р - давление жидкости; p - плотность;
F = (Fx,Fy,Fz)T
вектор массовых сил;
R = - вектор силы сопротивления, ко-
торое испытывают частицы жидкости в порах; п = (V хс, V ус, уи.)’г - вектор средней скорости частиц некоторого элементарного объема. Обозначив через V = (V х, V у, V ^г вектор скорости фильтрации, а через п - пористость грунта, можно написать:
V = ^ с. (2)
Считая, что ось Ое направлена вертикально вверх, имеем
? = (0, 0, - м)г . (3)
Использовав закон Дарси как линейный закон сопротивления, где последнее считается пропорциональным первой степени скорости фильтрации, П.Я. Полубаринова-Кочина [1] предложила зависимость
R. =-
gv x
Rv =-
gV V
R, =-
gv z
к у к ъ к ’ где к - коэффициент фильтрации. Тогда из уравнения (1) с учетом (2), (3) и (4) на основании вышеперечисленных допущений из гидродинамического уравнения Эйлера (1) получим линеаризованное и осредненное уравнение вида:
1 dv 1 , - g -
-------= — gradp + F + —v
n dt p k
(5)
П.Я. Полубариновой-Кочиной показано, что при значении коэффициента фильтрации k = 1 +100 (м/сут) членами, содержащими дифференцирование по времени, можно пренебречь, и тогда из уравнения (5) вытекает, что
p
(б)
V = -kgradH, Н = — + ъ.
РЕ
Здесь Н - пьезометрический напор; z - геометрическая высота. Следовательно, при перечисленных допущениях из уравнения Эйлера вытекает линейный закон Дарси (6).
Уравнения движения Эйлера содержат четыре неизвестные функции: Vх, Vу , у ъ и р. Присоединяя к этим уравнениям уравнения неразрывности с учетом несжимаемости жидкости и скелета грунта:
dvx dv y dv z
= 0,
(7)
Эх Эу Эъ
приходим к замкнутой системе (6), (7). Подставляя в (7) выражения для составляющих скорости фильтрации согласно (6)
, 9H
v x =-k ax'
, dH
v v =-k aT'
, dH
v z=-k sT
приходим к уравнению:
a
dx
aH
dx
_a_
Эу
k
dH
aV
dz
dH
k
dz
= 0 .
(8)
(9)
Если в пласте существуют источники и стоки подземных вод и возможен переток из выше-и нижележащих водоносных горизонтов, то
d
+
поступление из них в рассматриваемый пласт равно величине [2]
к к
8 = -^-(Ив - Н) + -^(Нн - Н) + f = ОН + f , тв тн
где: О = -
кв
тв
кн
тн
кк
W =^-Ив + ^Ин + f , т т
Нв, Нн - напоры в выше- и нижележащих напорных пластах;
к в, к н, и шв, тн - коэффициенты фильтрации и мощности верхней и нижней слабопроницаемых перемычек;
Г - функция источников и стоков. Прибавляя величину 8 в левую часть основного уравнения (9) и перенеся величину W в правую часть, для установившегося режима фильтрации подземных вод получаем уравнение:
Э С, ЭИ 4 к—
Эх
Эх
Э
ЭУ
к
ЭИ
ЭУ
Э
Эz
к
ЭИ
Эz
+ ОН = W,
(x,y,z) є V . (10)
Допустим, что V - цилиндрическая область, ограниченная поверхностью ЭV= Е = Ея + Е + Е ,
* * бок в и’
где Ебок - боковая поверхность цилиндра, Ев и Ен -верхнее и нижнее основания. Краевые условия на границе области имеют вид:
Эн
к ЭП - в1И = а1, (х, у, 2)єЕ(
бок,
к
ЭИ
(11)
(12)
(13) (13) в
(14)
Эп = 0 (х, у, 2)еЕн,
Н = Н1(х, у, 2), (х, у, 2) еЕв. Запишем граничные условия (11) общем виде:
У:к ^ЭН -у 2рН = У з«,
Эп
где:
У1 = 1> У2 = 1> Уз = 1> а = а1, Р = Р: на
У1 = 1, У2 = °, У3 = °, а = а1 на Ен,
У1 = 0, У2 = 1, Уз = 1, а = Н , р = 1 на Ев, где Ь = Ь (х, у, 2), а =а (х, у, 2), Н1(х, у, 2) - заданные функции.
Область V разбиваем на т треугольных прямых призм, тогда поверхности Ев и Ен также разобьются на т треугольников. Обозначим призму (фрагмент) через £ и приближенное решение задачи (10) - (13) ищем в виде:
ш п
Нп(х,у,ъ) = ЦН^х^ъ) = IН^(х,у,ъ), (15)
f=1 (Г) 1=1
где 1 = (х.,у.,2.) - глобальные узлы области V (количество их равно п); Р.(х, у, 2) - линейные базисные функции.
Применяя к задаче (10) - (14) обобщенный принцип Галеркина, имеем равенство [3]
№
Эх
\
+ОИп-W
Э С к ЭИп 4
Эх
Э
ЭУ
к ЭИп
ЭУ
Э Ск ЭИп 4
Эz
Эz
МШШк -ьи-“^“=0
Эп
І = 1,2,...,п, (16)
откуда, согласно формуле Грина, получаем:
ш
V
ЭК ЭН ЭКЭН ЭКЭИ
Эх Э :
п -
Э у Э у Э z Э ;
+№
ёу-
-я <ь+Яр)к1НПкйу -
-((Я,вН „ау -яявс1у =0, ]=1, 2,..., п. (17)
I I
Подставляя вместо функции Нп ее разложение (15) в равенство (17), имеем систему линейных алгебраических уравнений
(18)
где
ап = а„=Шк
'3^ +ЭЕ,
Эх Эх Эу Эу Эz Эz
ІУ +
+Ш3^ -ЦЯ^о, Ь =Я^\Щу+ЯЯ ado. (19)
V £ V X
Переходим теперь к описанию базисных функций.
Рассмотрим фрагмент, образованный из треугольной призмы с вершинами 1, _), к, 1' = 1 + А 1, / =) + А к' = к + Ак (рис. 1). Каждая такая элементарная призма состоит из трех тетраэдров: укк', к'/И и у ^к'.
к
Рисунок 1.
+
+
1=1
V
/
/
(2O)
В каждом тетраэдральном элементе (e) функция H(e)(x, y, z) образуется в виде
H(e) (x, y, z) = H(e)F (x, y, z) + H(e)F (x, y, z) +
+ Hke)Fk(x,y,z) + H(e)Fj(x,y,z), где i, j, k, l - вершины тетраэдра,
F (x,y,z) = (ai + bix + c;y+
Fj (x,y,z) = (aj + bjx+j + djzX
Fk(x,y,z)=(ak + bkx+cky+dkz),
Fj(x,y,z) = (at + btx+cty+djz),
C=
1x
1x
Vi Zi
i Vi Zi
1 xk Vk zk
1 х1 у1 ъ1 ,
Коэффициенты функций БДх, у, /)(8 = 1, .ь к, 1) являются элементами матрицы, обратной по отношению к матрице С, т. е. они являются алгебраическими дополнениями соответствующих элементов матрицы СТ, транспонированной по отношению к матрице С:
C-1 =
detC
ai ai bi ь. cc
I i dd
bk b1
dk d1
Н(е) = ^(л^у^Л § = 1,j,k,1.
Непосредственной проверкой можно установить, что функции формы Б;, FJ Бк и F1 равны единице в соответствующих вершинах 1, .), к, 1 и нулю - в других вершинах и что в произвольной точке элемента сумма их равна единице:
F1 (х, у,ъ) + К|(х,у,ъ) + Fk(x, у,ъ) + F1(x,y,z) = 1.
Поскольку область V аппроксимируется совокупностью 3т тетраэдральных элементов, то приближенное решение краевой задачи (10) - (13) ищется в виде
Нп(х,у,ъ) = 1 IН(е)(х,у,ъ) = £ [ H1(e)Fl(x,y,z)+
е = 1 (е) е = 1
+ Н(е)Е (х, у, ъ) + Н^ (х, у, ъ) + Н(е)Ц (х, у,ъ)].
Теперь выпишем явный вид коэффициентов и правых частей системы (18). Для каждого элемента (тетраэдра) с вершинами 1, .), к, 1 вычисляются шесть коэффициентов (по числу ребер):
а^ = k(e)(blbj + + dldj) + Q(e)-Р(Ач),
а1к = к(е)(Ь1Ьк + е1ек + dldk)+Q(e)-Р(£1к),
а11 = к(е)(Ь1Ь1+е1е1 + dld1)+Q(e) -р(А11), ajk = k(e)(bjbk + е^к + djdk)+Q(e)-р(А jk),
aj1 = к^^Ц+еjе1 + djd1)+Q(e)-р(А jl), ак1 =к(е)фЛ+екС1 +dkdl)+Q(e)-Р(А к1), где: к(е) = кС) • |аегС|/6,
Q(e) = ^ С/12°, р(А) = РС£) • 8а /12 ,
кС^ и Q(e|^ - средние значения соответственно функций к(х, у, 2) и Q(x, у, 2) в элементе (е). В случае равенства индексов получаются диагональные элементы системы. Для правых частей системы получаем выражения:
Здесь
(e) =
b. = W(e)+a(А l).
з,
W(e) = W(;) • ^еЮ| / 24, а(£) = а(£) • БА /
Wc(p) - среднее значение функции W(x, у, 2) в элементе (е); аС£) и РС£) - средние значения функций а(х, у, 2) и Р(х, у,/) соответственно на грани (А) тетраэдра (е) (для граней, не являющихся граничными, функции а и Р равны нулю); 8А -площадь грани (А).
Для диагональных элементов получаем формулу:
(А.)
где
a.. = k(e)(b2 + c2 + d2) + Q(e)-R. L
ii v i i У ^i Hi
О(е) = рс^. | аеіе | /60, в(Д^ = вт БД/б.
Матрица системы (18) является симметричной с диагональным преобладанием, поэтому она легко решается одним из точных или приближенных методов.
Алгоритм решения задачи (10) - (13) отлаживался на следующем тестовом примере.
В цилиндрической области:
V = {х2 + у2 <0,25, 0^<0,4
ищется функция Н(х, у, 2) при следующих данных: к(х,у,ъ) = х + у + ъ + 2, Q(x,y,z)=2,
W(x,y,z) = 2(х2 + у2 + ъ2 +1) - 8(х + у + ъ + 2) + 4.
На поверхности цилиндра задается краевое условие:
к—=РН+а,
Эп
где на боковой поверхности
а(х,у,ъ) = 0, Р(х,у,ъ) = к(х,у,ъ) ;
1,25+ъ2
a
a
k
1
c
c
k
на верхнем основании цилиндра: а(х,у,г) = 0, Р(х,у^) = 0 ; на нижнем основании:
а(х,у,г) = 0, Р(х,у,г) = ~2—хут^- .
х +у +z +1
Основание цилиндра разбито на 56 треугольных элементов. Шаги по х и у меняются от 0,1 до 0,2. Шаг по ъ равняется А ъ = 0,2, т. е. цилиндрическая область заменяется областью, состоящей из 112 треугольных прямых призм высотой Ь = 0,2 или из 336 тетраэдров. Общее число узлов области (вершин тетраэдров) равно 111.
Точным решением задачи является функция
Н(х,у^) = х2 + у2 + z2 +1-.
В табл. 1 приведены точные и приближенные значения этой функции в узлах второй четверти круга
х2 + у2 < 0,25 (х < 0, у > 0),
Таблица 1. Точные и приближенные значения напоров
№ азлов 7=0
2 4 5 9 10 11 16 18 19
Точн. значен. 1,252 1,245 1,163 1,252 1,163 1,080 1,250 1,040 1,000
Приал. значен. 1,274 1,273 1,176 1,266 1,197 1,077 1,272 1,015 0,903
№ азлов 7=0,2
39 41 42 46 47 48 53 55 56
Точн. значен. 1,292 1,285 1,202 1,292 1,202 1,120 1,290 1,080 1,040
Приал. значен. 1,242 1,314 1,174 1,360 1,244 1,095 1,408 1,051 0,902
№ азлов 7=0,4
76 78 79 83 84 85 90 92 93
Точн. значен. 1,412 1,405 1,322 1,412 1,322 1,240 1,410 1,200 1,160
Приал. значен. 1,321 1,445 1,297 1,562 1,427 1,207 1,551 1,147 0,935
причем узлы 1, 2, 4, 9, 16, 38, 39, 41, 46, 53, 75, 76, 78, 83, 90 лежат на боковой поверхности цилиндра, а узлы 19, 56 и 93 - на его оси. Наибольшие отклонения от точного значения наблюдаются именно в этих узлах, так как треугольные элементы, прилежащие к оси цилиндра, имеют наибольшие размеры: А х = А у = А ъ = 0,2, а ребра тетраэдра являются диагоналями этих треугольников.
Список использованной литературы:
1. Полубаринова - Кочина П.Я. Теория движения грунтовых вод. - М.: Наука, 1977. - 664 с.
2. Полубаринова-Кочина П.Я., Пряжинская В.Г., Эмих В.Н. Математические методы в вопросах орошения. - М.: Наука, 1969. - 414 с.
3. Джаныбеков Ч. Моделирование гидрогеодинамических процессов с применением ЭВМ.: Фрунзе: Илим, 1989. - 183 с.
