CC BY
28
УДК 517.9+519.633
МАТЕМАТИЧЕСКОЕ И КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ НЕКОТОРЫХ ФИЛЬТРАЦИОННЫХ ТЕЧЕНИЙ
СИДОРОВ М.В., СТОРОЖЕНКО А.В.
Рассматривается задача расчета фильтрационного течения жидкости под некоторым гидротехническим сооружением. Предлагается приближенный метод ее решения, основанный на совместном применении метода R-функций и вариационных методов.
1. Актуальность исследования
Задачи расчета движения жидкостей в пористых средах имеют большое значение при проектировании различных гидротехнических сооружений. Для решения этого класса задач математической физики используются различные точные и приближенные методы: метод Фурье, вариационные методы, методы теории функций комплексного переменного, метод фиктивных областей. Однако для каждого из этих методов имеются существенные ограничения, связанные с геометрией области, поведением коэффициента фильтрации и прочее [1, 2]. В связи с этим разработка приближенных методов решения этого класса задач является актуальной.
2. Постановка задачи
Целью данной работы является построение эффективных численных методов расчета фильтрации под гидротехническими сооружениями. Рассмотрим задачу движения несжимаемой жидкости под гидротехническим сооружением (плотиной). На рис. 1 приведена условная, но достаточно общая, схема фильтрации. Здесь D — область фильтрации, Dj — подводная часть плотины (флютбет), dQ 4 — граница водонепроницаемой области, dQ 5 — шпунт.
Стационарная фильтрация несжимаемой жидкости описывается в рамках линейного закона Дарси уравнениями
divu = 0, (1)
u = -k(x)Vh , x є d , (2)
где u — скорость фильтрации; k(x,y) — коэффициент фильтрации; h = у +— — пьезометрический
напор; p — давление; р — плотность жидкости, а g — ускорение силы тяжести.
В плоском случае x = (x,y) из (1), (2) следует уравнение для напора:
± k(x)®
dx { 5x
су
k(x)
Sh
dy
= 0,
D
(3)
При численном исследовании задачи удобнее от уравнения (3) для напора h(x) перейти к уравнению для функции тока y(x, у) с помощью соотношений
ux
uy
k( ) ^h
-k(y у)—
ox
k(x, у)
ch
dy
dy , dx
Тогда для функции тока получим уравнение
(4)
(5)
dx
1
1
k(x,y) dx ) dy \ k(x,y) dy
= 0, (x,y) є D .(6)
Дополним его соответствующими граничными условиями. На проницаемых участках границы SQj и dQ 3 расчётной области D ставятся однородные условия второго рода:
* = 0. (7)
dn SQ1U9Q3
Условие (7), как следует из (4) и (5), соответствует постоянству напора на SQj и Ш3 . Водоупор Ш4 , флютбет Ш3 со шпунтом Ш5 водонепроницаемы,
поэтому нормальная составляющая скорости и на этих участках границы равна нулю. Это приводит к тому, что на Ш3, Ш4 и SQ5 функция тока принимает постоянные значения. Для определённости положим, что
Мзп4 = Q , (8)
Ю3 1Ш5 = 0 . (9)
Величина Q задаёт общий расход жидкости.
Основные трудности, возникающие при решении задачи (6)-(9), следующие: неоднородность грунта (коэффициент фильтрации k(x,y) зависит от координат), границы dQi, i = 1,...,4 — криволинейны, наличие шпунта SQ5, который может быть не единственен.
3. Метод решения
Для решения поставленной задачи (6)-(9) будем применять метод R -функций (пусть шпунт dQ 5 отсутствует). Выделим область Q (пересечение |x|< L c D), вне которой фильтрация незначительна. На боковых границах Q поставим условие
Мx=+L = 0 . (10)
Предположим, что известны функции Ю1 (x, y), ю2 (x, y) , ®3(x,y) , ю4 (x, y) , такие, что:
РИ, 2004, № 4
1) Oi(x,y)|dQi = 0, і = ;
2) ®i(x,y) > 0
(11)
3)
dn
при (x,y) efiu(UdQj)U{x = ±L}; (12)
j=i
6toi = і і=i,.'*4. (13)
SQi
Такие функции всегда можно построить с помощью аппарата теории R-функций [3].
Тогда функция Ф1 (x, у) =
Q®2 (x, у)
D12^=xx (юі®3) xx+ХУ ^юію^ Ху •
ox ox Оу Оу
Недостатком структуры (20) является её громоздкость и наличие двух неопределённых компонент Фі и Ф2 • Решение задачи (16)-(18) можно упростить, если учесть, что условие (18) — естественное.
Легко проверить, что функция
Ф2 (x, у) = Ф1^,у) - (21)
1 2 2
®1 -®2 -®3 -®4 л« —(L -x ) (2)
2L----------Dj2) Фі (x,y)
1
22
®2 -®4 ла---(L - x ) + Юі-Юз
2L
удовлетворят условиям (4) и (5). Тогда в задаче (3)-(5) сделаем замену:
y(x, у) = Ф 2 (x, у) + u(x,y), (22)
где u(x,y) — новая неизвестная функция. Это приводит к задаче
Ф'
(x, у) + ю4 (x, у) ла (L2 - x2)(14)
2L
обладает свойствами: фі „ = 0 ,
І ґ\ 5^2
фрSq4y{x=±L = Q , на 5Q1 udO3 принимает, вообще говоря, произвольные значения.
Значит, задачу (6)-(10) можно записать в виде
_д_ dx
(15)
Ню2 l№4 U{x=+L = Ф^1 Э02 U9Q 4 U{x=+L (16)
(17)
( 1 д ( 1 ау^ д Г 1 au ^ а Г 1 au ^
ч k(x,y) дх у н ау v k(x,y) су 2 = 0, V(x, у) еО , 3x ^ k(x,y) dx J ду ^ k(x,y) ду J
ай
*2 *4 ^
= о
=f
ulЮ2 U9Q4 U{x=+L} 0 :
9Q1U3Q3
В [3] показано, что краевым условиям
I - J f^u
uSQ1 I 5n
удовлетворяет структура
5u
<5n
= 0
где обозначено
au
ю2 (18) f(x,y) = -
d
9Q1U3Q3
5ф2 ^ , d j 1 5ф2
V(x,y) ей,
(23)
(24)
(25)
(26)
u = ю^1 H--1—— [у + Ф2ф2 _D(2)(ra^1)-D(2)ф-
Ф1 +Ф2
-а^Ф1 -аф]+ф, (19)
где Ф1, Ф2 — неопределённые компоненты, d(2) 5ф2 3 , 5й2 d
Dі — * т •
1 8x dx ду ду .
12 2
Заменяя в (19) ю1 на ю2 'ю4 Ла-(L -x ), ю2
2L
на ®1 -Ю3 , ф на ф1, при a = 0 , ф = 0 получаем структуру решения задач (15)-(17) в виде
dx ^k(x,y) dx J ду ^k(x,y) ду
Введём в рассмотрение оператор краевой задачи (21)-(24), действующий в L2(q) по правилу:
. д
Au =------
1 5u
д
1 5u
dx ^ k(x,y) dx J ду ^ k(x,y) ду с областью определения
DA={u|u є c^q)n c2 (A, u so = 0, fj-
(27)
5Q1 U5Q3
= 0}
(28)
^ = Ю 2 ' ® 4 Л
1 (L2 - x 2 )Ф 1 +
2 L
DT2 2 4
®1 • ®2 '®3 'ю4 ла 2L(L _ x )
+------------------2L---------x
1 2 2
ю2 'ю4 ла 2L(L _ x ) +Ю1 Ю3 x [Ф1 -Ф3 -Ф 2 - D(2)(ra2 -®4 ла 2-(L2 _ x2)^)
- d(2)Ф1І + Ф1,
где Фь Ф2 — неопределённые компоненты, РИ, 2004, № 4
(20)
где дО. = дО 2 U 90 4 U {x=+А .
Пусть оператор (27) — эллиптический невырожденный. В этом случае можно доказать [4], что оператор А—положительно определён, условие (24) — главное, (25) — естественное. Значит, функция u может быть найдена как минимум функционала (Au,u) - 2(u,f) .
Тогда структуру решения краевой задачи (23)-(24) можно взять в виде:
1 2 2
u = ю2-®4 Л« —(L -x )ф , (29)
где ф — неопределённая компонента.
59
1
Для аппроксимации компоненты ф воспользуемся методом Ритца. Замкнув множество d(a) в
метрике, порожденной нормой
(' - \2
=я
1
HA Q k(x,y)
5u
dx
^ du ^
2 ^
Ч^У J
dxdy
(30)
/
получим энергетическое пространство HA со скалярным произведением
[u,v] = (Au,v) = Я
1 ( du dv Su dv
\
q k(x, y) ^Sx dx dy dy
dxdy.(31)
u
Выберем систему координатных функций фі , ф2 , ... , фи, удовлетворяющую условиям:
1) ~ф; є HA , i = 1,...,n ;
2) при любом n система фі, ф2,..., фn линейнонезависима;
3) система ~ф1, (вф2,..., Sфп полна в HA .
Неопределённую компоненту ф представим в виде
п
Ф = Е ciФі . (32)
i=1
Согласно методу Ритца для коэффициентов с1, ..., сп , получим систему линейных алгебраических уравнений:
Zl®фі’5фjlci = fфj), і = 1,...,n, (33)
і=і
где обозначено 5(x,y) = ю2
1 ,т:
ю4 Ла 2^ (L
x2).
Если область q может быть (при подходящем выборе системы координат) заключена в горизонтальную полосу 0 < y < H , то в качестве координатной системы можно выбрать
Фmn(x,y) = Pm [L^Pn - 1j , (35)
где Pk(t) — полиномы Лежацдра. Они образуют полную ортогональную на [-1,1] систему
Pk(t)
1 dk 2k k! dxk
[(x2
1)k].
(35)
N
Тогда u(x,y) = E cmnmn(x,y) . (36)
m+n=0
4. Результаты вычислительного эксперимента
Вычислительный эксперимент был проведён для случаев, когда шпунты отсутствуют, а флютбет имеет прямоугольное сечение. Для расчётов выбрано L = 3 , а высота водопроницаемого слоя H = 1. Рассмотрены случаи различного заглубления флют-бета (h = 0,25 и h = 0,5) с постоянной шириной, равной 1. Было принято, что границы 8&1 и SQ3 — плоские. Рассмотрены случаи, когда водоупор 3Q4 — плоский и неплоский. Во всех случаях коэффициент фильтрации принимался постоянным k(x, y) = 1, Q = 1.
1 случай. Заглубление флютбета h = 0,25 , водоупор плоский. В этом случае
ю1 -ю3 = 1 - y , ю2 = -(у - 0,75)ла (0,25 - x2),
Ю4 = y , d(2) =~ ,
dy
) _(у_0,75)ла (0,25-х2)
T1(x,y) =-------------------------------і-------;
- (у - 0,75) ла (0,25 - х2)+у ла 9 (9 - х2)
92(x,y) = Фl(x,y) +
+ (1 - уНу - 0,75) ла (0,25 - ^ ).У ла 9(9 - ^^
(у - 0,75) ла (0,25 - х2) • у ла 9(9 - х2) + (1 - у) ду
На рис. 2 приведена картина линий уровня (5ф = 0.1) приближённого решения для N = 5 (21 базисная функция).
і
0.8 0.6 0.4 0.2 0
Рис. 2. Линии уровня функции тока
-3-2-10 1 2 3
2 случай. Заглубление флютбета h = 0,5, водоупор плоский. В этом случае
Ю1 -®3 = 1 _y, ® 2 = Чу _0,5)ла (0,25 - x2),
Ю4 = y Ф1(x,y)=
d(2) =-—
д
dy
- (у - 0,5) ла (0,25 - х2)
- (у -0,5) ла (0,25 - х2) + у ла 9(9 - х2) ;
Ф2 (x,y) = фД^у) +
(1 - у) • (у - 0,5) ла (0,25 - х2) • у ла 9 (9 - х2) ^
(у - 0,5) ла (0,25 - х2) • у ла 9(9 _ х2) + (1 - у) ду
На рис. 3 приведена картина линий уровня (5ф = 0,1) приближённого решения для N = 5 (21 базисная функция).
0 L--- -----------l.
-3-2-10 1 2 3
Рис. 3. Линии уровня функции тока
3 случай. Заглубление флютбета h = 0,25 , водоупор неплоский. В этом случае
ю1 -ю3 = 1 - У , ю2 = -(у - 0,75Ла (0,25 - x2),
1
Ю4 = y-
Фl(x,y)=
= У--------x------ТЭ(2) —
1
-x 24 8
D(2) = —
’ ^1
dy
- (у - 0,75) ла (0,25 - х2)
- (у - 0,75) ла (0,25 - х2) + | у - -4-x - 8 1ла 1(9 - х2)
60
РИ, 2004, № 4
Ф2<Ху) = Фі(х,у) +
(1 - у) • (у - 0,75) л, (0,25 - х2) -(у - 2-х -8]ла 1(9 - х2)
(У - 0,75) л„ (0,25 - х2) у - 2-х - 8 8(9 - х2)+(1 - у)
На рис. 4 приведена картина линий уровня (5ф = 0.1) приближённого решения для N = 5 (21 базисная функция).
-3-2-10123
Рис. 4. Линии уровня функции тока
5. Выводы
Сравнение с аналогами. Полученные приближенные решения сравнивались с решениями, полученными П.Н. Вабищевичем в [1, 5] с помощью метода фиктивных областей. Результаты согласуются.
Преимущества предлагаемого метода по сравнению с существующими следующие:
— более высокая точность по сравнению с методом фиктивных областей (поскольку геометрия области учитывается точно);
— решение получается в аналитическом виде. Это облегчает его дальнейшее использование при на-9ф1 хождении поля скоростей и различных техничес-■ ких характеристик гидросооружения, чем определяется практическая значимость исследования.
Литература: 1. Лаврентьев М.А., Шабат Б.В. Методы теории функций комплексного переменного. М.: Наука, 1973. 736 с. 2. Вабищевич П.Н. Метод фиктивных областей в математической физике. М.: Изд-во МГУ, 1991. 156 с. 3. Рвачев В.Л. Теория R -функций и некоторые ее приложения. К.: Наук. думка, 1982. 552 с. 4 Михлин С.Г. Вариационные методы в математической физике. М.: Наука, 1970. 511 с. 5. Вабищевич П.Н, Гассиев Р.В. Численное решение задач напорной фильтрации под гидротехническим сооружением // ЖВМ и МФ, 1987. Т. 27, № 10. С. 1580 - 1584.
Поступила в редколлегию 21.11.2004
Рецензент: д-р физ.-мат. наук, проф. Дикарев В.А.
Сидоров Максим Викторович, ассистент кафедры прикладной математики ХНУРЭ. Научные интересы: математическое моделирование, математическая физика, теория R-функций и ее приложения. Увлечения и хобби: история культуры. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (0572) 702-14-36.
Стороженко Александра Владимировна, канд. техн. наук, доцент кафедры экономической кибернетики ХНУРЭ. Научные интересы: математическое моделирование экономических и технических систем, стохастический анализ и его приложения. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (057) 702-14-36.
РИ, 2004, № 4
61
