В.А. Хямяляйнен Установившаяся двумерная фильтрация жидкости вокруг перемычки

10

Та6лица 1

	Номер	PN	RN	Rb	h0b	hT	hy	lL	lT	lP	l	LL	LP	k0	k0b	КТ
	варианта	МПа						м							10-12 м2
	1	5	15	3	0,1	7	1	-	5	0,5	3	10	20	10	0,1	0,1
	2	5	15	3	0,1	7	1	-	5	1	3	11	19	20	0,09	0,09
	3	5	15	3	0,1	7	1	-	5	1,5	3	12	18	30	0,08	0,08
	4	5	15	3	0,1	7	1	-	5	2	3	13	17	40	0,07	0,07
	5	5	15	3	0,1	7	1	-	5	2,5	3	14	16	50	0,06	0,06
	6	5	15	3	0,1	7	1	-	5	3	3	15	15	60	0,05	0,05
	7	5	15	3	0,1	7	1	-	5	3,5	3	16	14	70	0,04	0,04
	8	5	15	3	0,1	7	1	-	5	4	3	17	13	80	0,03	0,03
	9	5	15	3	0,1	7	1	-	5	4,5	3	18	12	90	0,02	0,02
	10	5	15	3	0,1	7	1	-	5	5	3	19	11	100	0,01	0,01
	11	10	10	2,5	0,2	5	2	0,5	6	0,5	5	5	15	10	0,1	0,1
	12	10	10	2,5	0,2	5	2	1	6	1	5	6	14	20	0,09	0,09
	13	10	10	2,5	0,2	5	2	1,5	6	1,5	5	7	13	30	0,08	0,08
	14	10	10	2,5	0,2	5	2	2	6	2	5	8	12	40	0,07	0,07	10
	15	10	10	2,5	0,2	5	2	2,5	6	2,5	5	9	11	50	0,06	0,06
	16	10	10	2,5	0,2	5	2	3	6	3	5	10	10	60	0,05	0,05
	17	10	10	2,5	0,2	5	2	3,5	6	3,5	5	11	9	70	0,04	0,04
	18	10	10	2,5	0,2	5	2	4	6	4	5	12	8	80	0,03	0,03
	19	10	10	2,5	0,2	5	2	4,5	6	4,5	5	13	7	90	0,02	0,02
	20	10	10	2,5	0,2	5	2	5	6	5	5	14	6	100	0,01	0,01
	21	15	8	2,2	0,3	5	3	1	5	1	7	15	15	10	0,1	0,1
	22	15	8	2,2	0,3	5	3	2	5	2	7	16	14	20	0,09	0,09
	23	15	8	2,2	0,3	5	3	3	5	3	7	17	13	30	0,08	0,08
	24	15	8	2,2	0,3	5	3	4	5	4	7	18	12	40	0,07	0,07
	25	15	8	2,2	0,3	5	3	5	5	5	7	19	11	50	0,06	0,06
	26	15	8	2,2	0,3	5	3	6	5	6	7	20	10	60	0,05	0,05
	27	15	8	2,2	0,3	5	3	7	5	7	7	19	9	70	0,04	0,04
	28	15	8	2,2	0,3	5	3	8	5	8	7	18	8	80	0,03	0,03
	29	15	8	2,2	0,3	5	3	9	5	9	7	17	7	90	0,02	0,02
	30	15	8	2,2	0,3	5	3	10	5	10	7	l6	6	100	0,01	0,01

11

ПРИМЕРЫ ВЫПОЛНЕНИЯ ЗАДАНИЯ

Пример 1. Расчет распределения давления на прямоугольной области фильтрации

Постановка задачи. Определить величину притока воды через прямоугольный пласт песчаника с проницаемостью k1 = 10 Дарси, высотой H1=10 м, шириной D1 = 200м, длиной L1 = 100м. В середине пласта находится включение с проницаемостью k2 = 1 Дарси, с геометрическими параметрами: ширина L2 = 20 м, высота H2 = 7м (рис. 5). Оттоком с передней и задней стенками, а также в верхние и нижние слои пренебречь. Внешнее давление воды P0 = 105 Па.

	P=P	P=0
	y
H1		D
H2
	L2	x
z	L1

Рис.5. Установившееся течение воды через пласт

Решение. Так как течение установившееся, то распределение давления будет удовлетворять ДУ вида

∂∂ k( x, y,z ) ∂∂P + x x

	∂		∂P			∂		∂P		= 0 ,	(9)
		k( x, y,z )			+		k( x, y,z )
			∂y					∂y
	∂y					∂y

где x, y, z – декартовы координаты (рис. 5), k(x, y, z) – коэффициент проницаемости, м2.

12

В соответствии с симметрией задачи изменением давления вдоль оси z можно пренебречь. Уравнение (10) примет вид

∂∂ k( x, y,z ) ∂∂P + x x

	∂		∂P		= 0 .	(10)
		k( x, y,z )
			∂y
	∂y

Граничные условия в соответствии с постановкой задачи выглядят следующим образом:

∂P = 0,	x [0, L ], y = 0,		y = H	1	, так как нет перетока в верх-
∂y		1
ний и нижний пласты;
P = P0 ,	x = 0,	y = [0, H1 ], давление на левой границе;
P = 0,	x = L1,	y = [0, H1 ], давление на правой границе.

Нахождение распределения давления. В среде MatLab набираем pdetool (Enter).

Задаем масштаб рабочей области в появившейся графической оболочке PDETOOLBOX, выбрав в меню Options→ Axes Limits. В появившемся диалоговом окне задаем X-axis range [0 100], Y-axis range [0 10], нажимаем Apply, Close.

Рисуем область фильтрации, выбрав Draw→Rectangle/Square, начиная с точки с координатами (0, 0) и заканчивая точкой (100, 10).

Рисуем область включения, нажав первую кнопку на панели инструментов (прямоугольник) и начав рисовать от точки (40, 0) до точки (60, 7). Рабочая область будет выглядеть как на рис. 6. Сохраним про-

ект, выбрав File →Save As.

Рис.6. Область фильтрации в PDETOOLBOX

13

Область фильтрации можно также нарисовать с помощью команд

MatLab:

>> pderect([-20 30 2.2 12]); >>pderect([-20 30 2.2 3]).

Переходим в режим задания коэффициентов ДУ. Для этого выби-

раем PDE→PDE Mode, PDE→Show Subdomain Lables.

Зададим значения коэффициентов ДУ для каждой области течения. Коэффициенты приведены в табл. 2.

			Таблица 2
Область	с	a		f
1	10 10-12	0		0
2	1 10-12	0		0

Зададим граничные условия. Для этого необходимо перейти в режим задания граничных условий нажатием кнопки 7 панели инструментов и командой меню Boundary→ Show Edge Lables. Изменяем масштаб изображения, задав auto в диалоге масштаба. Двойным щелчком мыши на каждой границе задаем вид граничного условия и коэффициенты. Соответствующие значения приведены в табл. 3.

					Таблица 3
Граница	Вид граничного	g	q	h		r
	условия
1	Neumann	0	0	-		-
2	Neumann	0	0	-		-
3	Dirichlet	-	-	1		1500000
4	Dirichlet	-	-	1		1500000

После задания граничных условий и коэффициентов ДУ запускается расчет сетки и решения на области фильтрации нажатием кнопки 9 панели инструментов. Рабочая область будет выглядеть как на рис. 6.

Изменим параметры графика, выбрав Plot→ Parameters в меню и отметив галочкой пункт Arrous и в ниспадающем списке напротив – c*grad(u). После нажатия кнопки Plot рабочая область примет вид, изображенный на рис. 7.

По завершении полученную программу необходимо сохранить, выбрав File →Save. Код полученной программы можно посмотреть, от-

14

крыв соответствующий файл. Повторный запуск программы можно произвести, набрав имя файла в рабочей области MatLab и нажав Enter. Листинг получившейся программы приведен в прил. 1.

Рис. 7. Распределение давления и поле скоростей течения

Расчет притока воды. После нахождения распределения давления на области течения найдем приток воды в массив на границе 1 и отток из него через границу 4 (рис. 4). Для этого произведем экспорт решения и параметров разбиения области фильтрации в рабочую область

MatLab с помощью команд Solve→ Export solution и Mesh→ Export Mesh. В появившихся окнах диалога принять предлагаемые параметры.

Приток воды рассчитываем по формуле

	H	k(x, y)	∂P dy ,
	Q = D ∫1 −			(11)
	0	µ ∂x

где Q – приток, м3/с.

Таким образом, для нахождения притока необходимо найти градиент давления и вычислить интеграл (11). Расчет градиента давления производится с помощью функции gradient(uxy,dx,dy), где uxy – решение, заданное на прямоугольной сетке, dx – шаг сетки вдоль оси x, dy – шаг сетки вдоль оси y. Для нахождения значений решения на прямоугольной сетке необходимо создать сетку и вычислить решение на ней с помощью функции перехода от треугольной сетки к прямоугольной – tri2grid(p,t,u,x,y), где p, t – параметры треугольной

15

ной – tri2grid(p,t,u,x,y), где p, t – параметры треугольной сетки, u – распределение давления на области фильтрации, x – вектор-столбец разбиений вдоль оси x, y – вектор–столбец разбиений вдоль оси y. Матрицы p, t, u были импортированы в рабочую область MatLab ранее, векторы x и y необходимо создать. С учетом изложенного команды для нахождения притока на левой границе будут выглядеть:

>>k=10*10^-12;% коэффициент проницаемости >>mv=1.17*10^-3; % коэффициент динамической вязкости >>dx=0.01; % разбиение по оси x

>>dy=0.01;% разбиение по оси y >>x=0:dx:dx*5; % массив точек разбиения >>y=0:dy:10;% массив точек разбиения вдоль y

>>uxy=tri2grid(p,t,u,x,y); % расчет значений решения на полученной новой сетке

>>[px1,py1]=gradient(uxy,dx,dy);% расчет градиента давления

>>px=px1(1,:); % выбираем первый столбец для градиента вдоль оси x

>>pxv=px*k/mv; % скорость вдоль оси x

>>Q=-D*sum(pxv*dx)

Для нахождения интеграла (10) была использована аппроксимация интеграла конечной суммой с помощью функции MatLab sum(v). Аналогично находится величина притока на правой границе.

Пример 2. Найти скорость распределения и приток в выработку для геометрии в области фильтрации, представленной на рис. 1,при следующих параметрах :

LL = -20; LP = 30; Rb = 2,5; hоб = 0,7; RN = 12; l = 5; lL = 2,5; lT = 5; lP = 2,5; hy = 2; hT = 4; PN = 15 MPA; a = 100; в = 0,5; коб = 0,01; кN = 0,01.

Область фильтрации можно также нарисовать с помощью нескольких команд:

>> pderect([-20 30 2.2 12]); >>pderect([-20 30 2.2 3]);

>>pdepoly([-5.5 -5.5 -2.5 -2.5 2.5 2.5 5.5 5.5],[2.2 5 5 7 7 5 5 2.2]); >>pdepoly([-2.5 -2.5 0 0 2.5 2.5],[2.2 3 3.5 3.2 3.5 2.2]).

Что означают эти команды и их параметры можно посмотреть, на-

брав help pderect и help pdepoly в среде MatLab.

1. Переходим в режим задания коэффициентов ДУ. Для этого вы-

бираем PDE→PDE Mode, PDE→Show Subdomain Lables. Рабочая об-

16

ласть будет выглядеть как на рис. 6.

2. Задаем значения коэффициентов ДУ для каждой области течения. Коэффициенты приведены в табл. 4.

			Таблица 4
Область	с	a		f
1	y 0,1 10-12	0		0
2	y 0,01 10-12	0		0
3	y 0,01 10-12	0		0
4	y 100 exp(-0.5 y) 0,01 10-12	0		0
5	y 0,01 10-12	0		0
6	y 0,01 10-12	0		0
7	y 0,01 10-12	0		0
8	y 0,01 10-12	0		0

3. Задаем граничные условия. Для этого нажимаем кнопку 7 панели инструментов и выбираем Boundary→ Show Edge Lables. Изменяем масштаб изображения, задав auto в диалоге масштаба. Рабочая область будет выглядеть как на рис. 7. Двойным щелчком мыши на каждой границе задаем вид граничного условия и коэффициенты. Соответствующие значения приведены в табл. 5.

					Таблица 5
Граница	Вид граничного	g	q	h		r
	условия
1	Neumann	0	0	-		-
17	Neumann	0	0	-		-
15	Dirichlet	-	-	1		1500000
19	Dirichlet	-	-	1		1500000
6	Dirichlet	-	-	1		1500000
7	Dirichlet	-	-	1		1500000

4. После задания граничных условий и коэффициентов ДУ запускается расчет сетки и решения на области фильтрации нажатием кнопки 9 панели инструментов. Рабочая область будет выглядеть как на рис. 8.

5. Изменим параметры графика, выбрав Plot→ Parameters в меню

17

и отметив галочкой пункт Arrous и в ниспадающем списке напротив – c*grad(u). После нажатия кнопки Plot рабочая область примет вид, изображенный на рис. 9.

По завершении работу необходимо сохранить, выбрав File →Save. Код полученной программы можно посмотреть, открыв полученный файл.

Полный листинг полученной программы приведен в приложении.

Расчет притока воды в выработку. Расчет притока воды в выра-

ботку производится с помощью специальных функций MatLab. Для их работы необходимо импортировать данные из оболочки PDETOOLBOX в среду MatLab следующими командами PDE →Export PDE Coefficientы, Solve→ Export Solution, Mesh→ Export Mesh. Во всех диа-

логовых окнах принимать предложенные имена переменных. Пример подпрограммы для вычисления притока воды в выработку на участке границы 15 и 16 представлен в листинге 1.

%///////////// Листинг 1.

function inFl=inflowFirstWIF(xmin,xmax,ymin,p,t,u,k)

%функция расчета притока для границы вдоль выработки

%в секунду на метр периметра выработки

%xmin – минимальная координата участка по оси x

%, xmax – максимальная координата участка по оси x

%ymin – радиус, на котором вычисляется водоприток %p,t – параметры разбиения области

%u – решение ДУ

%k – коэффициент проницаемости

%mv – коэффициент динамической вязкости воды

global mv

dx1=0.025; % разбиение по оси x dy1=0.025;% разбиение по оси y x=xmin:dx1:xmax; % массив точек разбиения

y=ymin:dy1:ymin+dy1*2;% массив точек разбиения вдоль y I=-(xmin-xmax)/dx1; % количество точек разбиения y1=ones(I+1,1)*dx1; % вспомогательный массив uxy=tri2grid(p,t,u,x,y); % расчет значений решения на полученной новой сетке

[px1,py1]=gradient(uxy,dx1,dy1);% расчет градиента давления px=px1(1,:); % выбираем первый столбец для градиента вдоль

18

оси x

py=py1(1,:); % выбираем первый столбец для градиента вдоль оси y

pxv=px*k/mv; % скорость вдоль оси x pyv=py*k/mv; % скорость вдоль оси y inFly=pyv*y1; % суммируем inFlx=pxv*y1;

inFl=[inFlx inFly]; % возвращаем результат в виде матрицы

////////// Конец листинга 1

При расчете притока на границах 8 и 9 необходимо учесть, что расчет проводится в цилиндрических координатах. Пример подпрограммы для вычисления притока воды в выработку на данной границе представлен в листинге 2.

%//////Листинг 2.

function inFl=inflowThirdWI(xmin,xmax,ymin,ymax,p,t,u,a,b)

%расчет водопритока через массив global mv

dx1=0.05;

dy1=0.1;

x=xmin:dx1:xmin+dx1*3;

y=ymin:dy1:ymax; I=(ymax-ymin)/dy1; y1=ones(1,I+1)*dy1;

kk=a*exp(b*y);

%расчет притока слева

uxy=tri2grid(p,t,u,x,y);

uxy(I+1,:)=uxy(I,:);

[px1,py1]=gradient(uxy,dx1,dy1);

px2=px1(:,1);

py2=py1(:,1);

px2=kk'.*px2/mv;

19

py2=kk'.*py2/mv;

inFlefty=y*py2*dy1;

inFleftx=y*px2*dy1;

%расчет притока справа x=xmax-dx1*3:dx1:xmax; %y=ymin:dy1:ymax; %I=(ymax-ymin)/dy1; %y1=ones(1,J+1)*dy1;

uxy=tri2grid(p,t,u,x,y);

uxy(I+1,:)=uxy(I,:);

[px1,py1]=gradient(uxy,dx1,dy1);

px2=px1(:,3).*kk'/mv; py2=py1(:,3).*kk'/mv; inFrighty=y*py2*dy1; inFrightx=y*px2*dy1;

MQright=[inFrighty inFrightx]; MQleft=[inFlefty inFleftx]; %Mqdif=MQin-MQout; %MQdif=[Mqdif;Mqdif/inFrightx*100] inFl=[MQleft;MQright];

% Конец листинга 2.

СПИСОК РЕКОМЕНДУЕМОЙ ЛИТЕРАТУРЫ

1.Сборник инструкций и других нормативных документов по технике безопасности для угольной промышленности. – М.: Недра, 1978. – 744 с.

2.Хямяляйнен В.А. Формирование цементационных завес вокруг капитальных горных выработок / В.А. Хямяляйнен, Ю.В. Бурков,

П.С. Сыркин. – М.: Недра, 1996. – 400с.

3.Заславский Ю.З. Новые виды крепи горных выработок / Ю.З. Заславский, Е.Б. Дружко. – М.: Недра, 1989. – 256с.