Современные методы математического моделирования при прогнозе геологических нарушений угольных пластов Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

---------------------------------------- © А. А. Глухов, А.В. Анциферов,

2005

УДК 622.031:51.001.57

А.А. Глухов, А.В. Анциферов

СОВРЕМЕННЫЕ МЕТОДЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ПРИ ПРОГНОЗЕ ГЕОЛОГИЧЕСКИХ НАРУШЕНИЙ УГОЛЬНЫХ ПЛАСТОВ

Семинар № 2

ТЪ настоящей работе рассмотрены конечно-разностные соотношения, позво-

-Ш-9 ляющие моделировать процесс распространения сейсмических колебаний в угленосной толще. Предложены способы учета напряженного состояния горного массива, пористости, трещиноватости угля и пород при математическом моделировании процесса распространения сейсмических колебаний в угленосной толще и учета поглощения сейсмических колебаний в среде при использовании МКР.

Новейшие методики исследования нарушений угольных пластов основаны на предварительном математическом моделировании процесса распространения сейсмических волн [1-5]. Модель для расчетов выбирается на основе априорных данных

о горно-геологических условиях залегания пласта и о нарушении. Итоги моделирования принимаются в расчет при выборе оптимальных условий проведения натурных наблюдений и методов обработки их результатов. Применительно к рассматриваемой проблеме наиболее распространенными методами математического моделирования являются конечно-разностные, реализованные в ряде практических разработок. Он обладает существенными преимуществами перед аналитическими методами, но не лишен своих недостатков.

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

Второй недостаток основан на том, что ранее использованные подходы используют расчетные уравнения Ламе для однородной изотропной среды. Хотя данное упрощение во многом оправдано, для ряда моделей (особенно для тех, где присутствует градиентное изменение физико-механических параметров угля) следует применить иной подход.

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

Как известно, уравнение движения упругой среды в общем виде можно записать

как

догк

Р = 1^, (1

дхк

д&гк

где ------ представляет собой силу внутренних напряжений, р - плотность, и г - ус-

дхк

корение, О ш - тензор напряжения, который в соответствии с законом Гука можно записать в виде:

О = К„А + Ц^ -1^,,,,), (2)

2

где игк - компоненты тензора смещения, р - плотность среды, К = Л +— и представля-

3

ет собой модуль всестороннего сжатия, X и ^ - коэффициенты Ламе.

°к = Ли11^к + 2иик

С учетом этого, выражение (1) можно записать в виде:

дЛ е _ ди . е дип дигк

Ри* = д и^гк + 2игк + Л$гк дГ + 2и~д,------------------------------------------------- . (3)

дхк дхк дхк дхк

Если в выражение справа от знака равенства опустить первые два слагаемых, то получится известное уравнение Ламе, записанное для изотропной однородной среды и реализованное в [3, 4, 5]. Первые два слагаемых необходимы при учете неоднородности среды.

Основной проблемой при использовании конечно-разностных методов на ЭВМ является тот факт, что границы решетки порождают «паразитическую» часть решения, эквивалентную отраженной от свободной границы волне. В шахтной сейсморазведке при моделировании распространения сейсмических колебаний в пласте эта проблема осложняется тем, что скорость распространения сдвиговых волн в угле может в 5-6 и более раз быть меньше, чем скорость распространения продольных волн в породах. В результате на полезную часть решения могут накладываться многократные отражения от любых границ решетки.

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

—ах t

и = ие х ,

-ау t V = \е у ,

—а t

^ = т^е а ,

где аг - коэффициенты затухания компонент смещения во времени.

Тогда выражение (3) можно записать в виде системы:

. 2 \ /. \ д ( ди дv дм) ( д2 и д2и д2и )

рри — аи + ам )=\Л + и)—I-----------1----1-I + и1 —о- +------------о- +-I +

У x x ’ У дх ^ дх ду дz) удх2 ду2 дг ) + дЛ ( ди + ^ + дм ) + 2 ди ди + ди ( ди + дv) + ди ( ди + дм

дx ^ дх ду дz) дх дх ду ^ ду дх) дz ^ дz дх

. 2\ чд(ди дv дм) (д 2v д 2v д 2v) (4)

р—а+а)л+и) + * )+^[зх2V+д? 1+

+ дЛ( ди + дv + дм) + ^ ди дv + ди I ди + дv) + ди ( дv + дм

ду ^ дх ду дz) ду ду дх ^ ду дх) дz ^ дz ду

л. . 2 \ / \ д (ди дv дм) (д2м д2м д2м

рм—ам+ам)=(Л + и)\ 1Т+1Т+1Г| + и| -гг+1тт+!гт| +

+ дЛ ( ди + дv + дм ) + 2 ди дм + ди ( ди + дм I + ди ( дv + дм

дz ^ дх ду дz) дz дz дх \дг дх) ду ^ дz ду / где u,v и м - компоненты смещений по х, у и г соответственно.

Моделируемая область представляется неравномерной решеткой из М*Ы*Ь элементов (по осям х, г и у соответственно) размерами

Ахт,п,, X Агт,п,, х Аут,п,, каждый (т = 1..М, п = 1..Ы, I = 1..Ь). Величины Ахщп,,

и Аут,„:1 представляют собой шаги дискретизации модельной решетки в пространстве, Дt - во времени (номера шаговр = 1..Р). Расчетные соотношения можно записать в виде:

ир+1 = а ,р—1 + Ъ ,р + с ,р + С ,р +

т,п ,1 х т,п,1 т,п,1 х т,п,1 т,п ,Г ^'xm+0.5,„lrm+1,n^т‘'хт—0.5,п,Гт—1,п,1 ^

+ Сх т,п+0.5,{ит,п+\,1 + Сх т,п—0.5,{ит,п—\,1 + Схт,п,1—0.5ит,п,1—1 + Схт,п,1+0.5ит,п,1+1 +

1 хг р л хг р л хг р л хг р

+ т,п+\,1мт,п+\,1 т,п—\,1мт,п—\,1 + т+1,п,/мт+1,п,/ т—\,п,1мт—\,п,1 +

+ Иху л?р — Иху л?р + Их л?р — Их л?р +

т,п,,+1 т,п,,+1 т,п,,—1 т,п,,—1 т+1,п,, т+1,п,, т—1,п,, т—1,п,,

+ {х2 м?р + {Х2 м?р + {Х2 м?р +

т+0.5,п+0.5,, т+1,п+1,, т—0.5,п—0.5,, т—1,п—1,, т—0.5,п+0.5,, т—1,п+1,,

+ {хг мр + {ху Л)р + {ху Л)р +

т+0.5,п—0.5,, т+1,п—1,, т+0.5,п,,+0.5 т+1,п,,+1 т—0.5,п,,—0.5 т—1,п,,—1

т+0.5,п—0.5,, т+1,п—1,, т+0.5,п,,+0.5 т+1,п,,+1 т—0.5,п,,—0.5 т—1,п,,—1

+ {ху vp + vp

т—0.5,п,,+0.5 т—1,п,,+0.5 т+0.5,п,,—0.5 т+1,п,,—1

(5)

;р+1 = а vp 1 + Ъ vp + с vp + с vp

т,п,, у т,п,, т,п,, у т,п,, т,п,, у т+0.5,п,, т+1,п,, у т-0.5,п,, т-

рррр

+ с ^ + с ^ + с ^ + с ^ +

у т,п+0.5,, т,п+1,, у т,п-0.5,, т,п-1,, у т,п,,-0.5 т,п,,-1 у т,п,,+0.5 т,п,,+1

+ И у2 мр — И у2 мр + И у мр — Иу2 мр +

т,п+1,, т,п+1,, т,п-1,, т,п-1,, т,п,,+1 т,п,,+1 т,п,,-1 т,п,,-1

+ Иух +1ит ,+1 — Иух} мр , 1 + Иух и , — Иух 1 ,ит 1 , +

,п,,+1 ,п,,+1 ,п,, 1 ,п,, 1 +1,п,, +1,п,, 1,п,, 1,п,,

+ / т,п+0.5,1+0.5мт,п+\,1+1 + /т,п-0.5,,—0.5 мР,п—1,,—1 + /т,п+0.5,1—0.5 ^т,п+1,1—1 +

+ {у2 мр + {ух ,р + {ух ,р +

Jm,п—0.5,,+0.5 т,п—1,1+1 J т+0.5,п,1+0.5 т+1,п,1+1 У т—0.5,п,,—0.5 т—\,п,1—1

т,п—0.5,,+0.5 т,п—1,1+1 л т+0.5,п,1+0.5 т+1,п,1+1 ^ т-0.5,п,1—0.5 т—1,п,1—1

ух ,р + {ух ,р

0.5,п,,+0.5 1,п,,+0.5 +0.5,п,, 0.5 +1,п,, 1

+ {ух ,р + {ух ,р

0.5,п,,+0.5 1,п,,+0.5 +0.5,п,, 0.5

(6) +

мт,п,1 аг т,п,1мт,п,1 + Ъг т,п,{мт,п,1 + сг т+0.5,п,/мш+1,п,/ + Сгт—0.5,п,/мш—1,п,1 +

+ сг ,п+0.5,,м р,п+1,, + сг ,п 0.5,,м р,п 1,, + сг ,п,, 0.5м р,п,, 1 + сг ,п,,+0.5м р,п,,+1 +

. згх , р 1 гх „ р . 1 гх „ р згх „ р .

+ И .1 и ,17 — И 1 17+ И . 1 . 1 1 — Ил 1Ы л 1 +

,п+1,, ,п+1,, ,п 1,, ,п 1,, +1,п,, +1,п,, 1,п,, 1,п,,

+ Игу л?р —Игу л?р + Игу л?р —Игу л?р + (7)

,п,,+1 ,п,,+1 ,п,, 1 ,п,, 1 ,п+1,, ,п+1,, ,п 1,, ,п 1,,

+ {гх ,р + {2Х ,р + {2х ,р +

+0.5,п+0.5,, +1,п+1,, 0.5,п 0.5,, 1,п 1,, 0.5,п+0.5,, 1,п+1,,

+ {гх и р + {г vp + {г vp +

+0.5,п 0.5,, +1,п 1,, ,п+0.5,,+0.5 +1,п,,+1 ,п 0.5,, 0.5 ,п 1,, 1

+ {гу vp + vp

,п 0.5,,+0.5 ,п 0.5,,+0.5 ,п+0.5,, 0.5 ,п+1,, 1

Полагая коэффициенты затухания в зонах у границ решетки одинаковыми для всех компонент (а=ах= ау= а^), что ни коим образом не влияет на решение задачи, имеем:

(

а , = а , = а ,

х т,п,1 у т,п,1 г т,п,1

1 +-

а

\

V 2А)

0 г! (Лт,п,1 + 2ит,п,1 )

ч

х т,п,1

= ч 2 — 2

Мр

Ъ

Ах2

п. п.( Л'т,п,1 + 2Мт,п,1 )

+ • т,п,1 + Мт ,п,1

т,п,1 Аг2 ,п,, Ау2 т, п,1 )

Р

т, п,1

— а

у т,п,1

= ч 2 — 2

п,1) + Мт, п,1 + Н1 т, п,1

\

(

Ау т,п,1

Му,

Аг 2 ,п,, Ах2 т,п,1 )

Р

т, п,1

а

=ч

2 — 2

(Лт ,п,1 + 2Мт ,п,1 )

п,1/ + Мт ,п,1 + и ,п,,

2

)

Ах2

РП.

а

Ъ

г т,п

,п,,

,п

,п

с = к

х т+0.5,п,, т,п,1 * 2

Ах т,п,1

Х (/! (Лт+1,п,, Лт- 1,п,1 )+ 2 )ит+1,п,1 Мт —1,п,1 )+ 2Мт+0.5,п,, + Л

+0.5,п,,

X

_ k 1

c с с j — —k і —----------------------x

xm—0.5,n,l m,n,l 2

AX m,n,l ,

X (4 (Лр+1,п,1 Лгп—l,n,l )+ 2 lUm+l,n,l Um—l,n,l ) 2—m—0.5,n,l Лр—0.5,n,l )

_ k 1

Czm,n+0.5,l m,n,l і 2 X

Az m,n,l ,

X )4 (,n+l,l — Лп,п—l,l )+ "2 )um,n+l,l — U m,n—l,l )+ 2—m,n+0.5,l + Лm,n+0.5,l )

_ k 1

c с с j — —k і — --------x

z m,n—0.5,l m,n,l 2

Az m,n,l ,

X ("4 (^'m,n+l,l Лш,п—l,l )+ 2 lUm,n+l,l Um,n—l,l ) 2—m,n—0.5,l Лш,п—0.5,l )

_ k 1

Cym,n,l+0.5 m,n,l д 2 X

Ay m,n,l ,

X (44 (m,n,l+1 — Лш,п,і—1 )+ "2 )-m,n,l+1 — -m,n,l—1 )+ 2—m,n,l+0.5 + Лп,п,і+0.5 )

_ k 1

c n с j — —k і — --------x

y m,n—0.5,l m,n,l 2

Ау m,n,l ,

X (4 (m,n,l+l Лп,п,1—1 )+ 2 )—m,n,l+1 Um,n,l—1 ) 2—m,n,l—0.5 Лш,п,1—0.5 )

_ k 1 U _ k 1 U

Cxm,n+0.5,l m,n,l д 2 -m,n+ 0.5,l , Cxm,n—0.5,l m,n,l д 2 —m,n-0.5,l ;

Az m,n,l Az m,n,l

_ k 1 ' _ k 1 ,

Cxm,n,l+0.5 m,n,l д 2 —m,n,l + 0.5 , Cxm,n,l—0.5 m,n,l д 2 -m,n,l—0.5 ;

Ay m,n,l Ay m,n,l

c k 1 U _ 1 1 '

Cy m+0.5,n,l _ km,n,l д 2 -m+ 0.5,n,l , Cym-0.5,n,l km,n,l Ax2 —m—0.5,n,l ,

Ax m,n,l Ax m,n,l

c k 1 . / _ k 1 '

Cym,n + 0.5,l _ km,n,l K 2 Um,n + 0.5,l, Cym,n—0.5,l m,n,l A 2 —m,n—0.5,l ,

Az m,n,l Az m,n,l

_ k 1 ' _ k 1 ,

Czm+0.5,n,l m,n,l д 2 + 0.5,n,l , Czm—0.5,n,l m,n,l д 2 Um—0.5,n,l ,

AX m,n,l AX m,n,l

_ k 1 ' _ k 1 ,

Czm,n,l+0.5 m,n,l д 2 Um,n,l+0.5, Cz m,n,l—0.5 m,n,l д 2 —m,n,-0.5l ,

Ay m,n,l Ay m,n,l

dxz _ 1 k ______________________1_________)// — / / ) dxz _ —dxz

m+1, n,l 4 m,n,l Ду Az \*m, n+l,l A^m, n—l,l/> m—l,n,l m+l,n,l J

m,n,l m,n,l

dxz 1 k 1 (л Л ) dxz dxz

m,n+l,l 4 m,n,l ж ж V^p+l,n,l ^m—l,n,l), m,n—l,l w,n+l,l -

Axm,n,lAzm,n,l

m,n,l m,n,l

dzx _ .1 k ____________________1__________)Л — Л ) )zx _—dzx

m+1,n,l 4 w,n,l ж ж V m,n+l,l m,n—1,1/, w—l,n,l w+l,n,l -

Ax l Az l

m,n,l m,n,l

fu‘m,ru‘m

i ^S'0+l‘“‘S'0-tUjj _|_ S0+r«‘S'0-“^_____________I___________£'o+f«‘s o-«y _ so+f«Vo-«y

AV -W < <

s^'O-/w _|_ &'0-/w 1&'0-/w &'0-«y _ S'O-/w &'0-«y

/‘m!w /„/‘m!w4„

AV -W < <

^&'0+/w &'0+«'^y _|_ £'0+/« £'0+«^_______________1__________Z W M/^±.— £'0+/« £'0+«y _ S‘0+/w &'0+«y

7‘m‘m/_i‘u‘ui _

Z\7 XV < < < <

i ^ I £'0~u £'0+mjj _|_ /S'O-w £'0+«^_________1___________y_lum^ P___ I S' 0-w &'0+«y _ /V0-« V0+«y

7‘h‘M/_/‘M!W ______

' zv xv

&'0+w &'0-«'^y _|_ / £'0+« £'0-«^_________1____________l u m^±._ — ?'£'0+« S'0-«y _ / S'O+w S'0-«y

7‘M‘M/__7‘h‘M/____

zv XV

&'0-w &'0-«'^y _|_ /S'O-w £'0-«^__________1___________W / &'0-w &'0-«y _ /Vo-w £'0-«y

7‘h‘M/__7‘h‘M/____

zv XV

&'0+« &'0+«'^y _|_ / £'0+« £'0+«^_________1___________W / £'0+« £'0+«y _ / V0+«V0+«y

7‘M‘M/ ./* 7‘U‘IU __

1 Ay1 zy

i+/

u‘mp_ — I~Vu‘mp ‘^lc\-u‘mrf _ /‘l+W MY/j___________________’___________l‘u‘m3jlL — l+l‘u‘mp

1 I + W _ 11" ""' » I I I '•* “^/ I 1 Z "-1 X \ ___ I "•' >/ k _ I I 1

1 I —

1+/ L --1_ _ L I -> — _!* <. I [ I ~ I L 1 " -1 X \ ________________ I ^ __ I 1 Z "•'

1 I —

/‘I+w

1+/

mp_ — l‘l~U‘mp ‘^\. VU‘myj _ l + l‘u ujff^_________________ _______________ j‘u‘Wyt_ _ fl+u‘mp

tUip_ — I ~l‘U‘mp ‘ _ VU‘\ + lUffj___________________ ____________ fu‘u4yb_ _ \ + fu‘U4p

1‘u‘i+tu^ _ 1‘u‘i-tu^ _ l+fu‘mv\________’__________y_l‘u‘WyV_ _ Vu‘l+mp

1+l‘u‘ta^ _ i6(l‘u‘l-^ _ ru‘l+lu1/\_________1_y_Vu‘w^\^ _ i+i‘u‘«ip

fu‘\+UA _ 1‘U‘I-UA I+/‘M‘K^y\__________;______________V f _ /‘W‘I+W

AxP — AxP ( j ^ 7 [ — AxP

fu‘iu___7‘M‘M/______

, l‘l+u‘iUp_ _ 1‘1-u‘iUp < ^ _ fu‘\+iUjY^__________V__________v_ fu‘iuyfr_ _ fl+u‘iUp

I I

/eM 1 m

1 I

/eM 1 "w,,

1 I

Zew I Ul ^ r.

1 I

/eW ‘“VCv'‘"‘ “w,,

1 I

y\ xy

I

xy

I

r«

1 I

ti 'Mzyl‘u w r

/ху

J т+0.5,п,/-0.5

/уг

J т,п+0.5,/+0.5

= /ух

J т+0.5,п,/-0.5

1

Ах.

= {гУ = 1 к

у т,п+0.5,/+0.5 4 т,п,/

т,п,/Аут,п,/ 1

(Л

Ау , Аг

у т, п,/ т

'т+ 0.5, п,/-0.5 + Мт+0.5, п,/-0.5, п+0.5,/+0.5 + Мт,п + 0.5,/+0.5 ),

,/

/

Уг

т,п-0.5,/-0.5 /т,п-0.5,/-0.5 4 ^т,п,1 ж ж \Лт,п-0.5,/-0.5 + Мт,п-0.5,/-0.5 /

Ау , А /

у т,п,/ т,п,/

руг

т

т,п+0.5,/-0.5

/т

уг

/т,п+0.5,/-0.5 4 ^т,п,1 А А \Лт,п+0.5,/-0.5 + Мт,п+ 0.5,/-0.5^

4 Ау , Аг,

у т,п,/ т,п,/

= -1 к _______________1________(л

4 т,п,/ ж ж \ т

т,п-0.5,/+0.5 /т, п-0.5,/+0.5 4 кт,п,/ * л \Лт,п-0.5,/+0.5 + Мт,п-0.5,/+0.5 >

Ау . Аг

у т,п,/ т

Мт+0.5,п,/ М Мт-0.5,п,/ М Мт ,п+0.5,/ М Мт ,п-0.5,/ = М Мт ,п,/+0.5 = М

т,п,/-0.5

= М

рт

м1

Рт,п

'т+0.5,п,/ + 4 (мт+1,п,/

т-0.5,п,/ -1 (м 4 т +1, п , /

т,п+0.5,/ + 4 (мт,п+1,/

т,п-0.5,/ -1 (М 4 V' т,п+1,/

т,п,/+0.5 + 4 (Мт,п,/+1

т,п,/-0.5 -1 (М 4 т,п,/+1

кт,п,/ т,п,/ Рт,п,гЧ , (

ч =

т,п,/ т,п,/

Мт-1,п,/

Мт-1,п,/

Мт,п-1,/

Мт,п-1,/

М т,п,/-1 М т,п,/-1

2А/

2 А -а

Выбор шагов дискретизации модели определяется характерными расстояниями, на которых изменяются характеристики угля и пород. Использование неравномерной решетки в данном случае является обоснованным, поскольку углевмещающая толща за исключением зон нарушений может быть представлена в виде слоисто-однородной среды. Скачкообразные либо градиентные изменения физических характеристик угля и пород в плоскостях, параллельных залеганию пласта наблюдаются только в зонах нарушения. Этот факт может быть учтен в геометрии решетки. Характеристики толщи в плоскостях, которые перпендикулярны залеганию пласта, следует рассматривать в зависимости от частотной характеристики поля колебаний как скачкообразно (от породы к породе) либо градиентно изменяющиеся во всем объеме модели. Алгоритм оптимального выбора топологии решетки таким образом, чтобы уменьшить общее число ее узлов (для уменьшения времени счета) и в то же время адекватно описать среду представляет собой отдельную задачу и не рассматривается в пределах данной статьи. Следует только отметить тот факт, что в настоящее время трехмерное моделирование требует большого объема вычислений, что является причиной его ограниченного применения.

Для плоской задачи, предполагая отсутствие зависимости по у, систему (4) можно записать в виде:

. 2 \ / л \ д ( ды ш

р\ы-ахы + ахы) =(л + м)— I ~дх + ~дг) + М

кдх2 дг2 у

+

+ дЛ ды + дм | + 2 дц ды + дц( ды + дw

дх I дх дг I дх дх дг I дг дх У

р\у -а,у\ + ауУ\ — ц

( д2у д2у^

чдх2 дг2 у

+

дц ду дц ду

дх дх дг дг '

?((

р\м> -а^ + а1м) —

)= (л + М)-^

д I ды дм> дг I дх дг

+ м

( д2^ д2^ ^

- + -

ч дх2 дг2 у

+

+ дЛ ( ды + дм ) + 2 дц дм + дц ( ды + дм

дг ^ дх дг) дг дг дх ^ дг дх

Моделируемая область, как и в трехмерном случае, представляется неравномерной решеткой из М хЫ элементов (по осям х и г соответственно) размером

Ах х Аг каждый (т=1..М, п=1..Ы). Расчетные соотношения можно записать в

т ,п т ,п ^ \ ’ /

виде:

ыр — а ыр + Ь ыр + с + п ^ ы'м + с г.с ы' 1 + с + п^ы' ,1 +

т,п хт,п т,п хт,п т,п хт+0.5,п т+1,п хт-0.5,п т-1,п хт,п+0.5 т,п+1

,Р

,р

+ с„ „-п Х,п-1 + ^ „+, wp „+, - ^ wp„-П + ^ ^р+,„ - ^ „-мр^„ +

хт,п-0.5 т,п-1 хт,п+1 т,п+1 хт,п-1 т,п-1 хт+1,п т+1,п хт-1,п т-1,п

т+0.5,п+0.5^т+1,п+1 + /т-0.5,п-0.5 ^т-1,п-1 + /т-0.5,п+0.5 ^т-1,п+1 + /т+0.5,п-0.5 ^т+1,п-\

(9)

^т,п + Ьгт,п^т,п + С г т+0.5,п^т+1,п + С г т-0.5,п^т-1,п + С г т,п+0.5^т,п+1 +

,Р

,р

+ с п ^м^ 1+ ^ ^ ыр .1 — ^ 1 ыр 1+ ^ ^ ы'. 1 — ^ 1 ы' 1 +

гт,п-0.5 т,п-1 гт,п+1 т,п+1 гт,п-1 т,п-1 гт+1,п т+1,п гт-1,п т-1,п

(10)

+ /

т+0.5,п+0.5 т+1,п+1

+ /т-0.5,п-0.5ыт-1,п-1 + /т-

ыР + / ыР

0.5,п+0.5мт-1,п+1 У т+0.5,п-0.5 т+1,п-1

уР+1 — а„,.„,ур- + Ь,.__уР_ + с

,р-1

р

+ с.

.у

р

у т,п т,п у т,п т,п у т+0.5,п т+1,п у т-0.5,п т-1,п

+

+с

5у

р

+ су

у т,п+0.5 т,п+1 у т,п-0.5 т,п-1

5у

р

где

а — а — а — -

хт,п ут,п г т,п

1 + -

а

2Ы

Ьх т,п — Ч

о 0((лт,п + 2цт,п )

2-2

Ах2

+

>

4. Ч , Ч

ГУ

Цт,п \

Аг т, пУ

2Ы

2А/ — а

\

Р<т,п -а ■■

У

(11)

Р

т,п

by m,n — q

2 - 2f ----- ----+------ ---- lp^,n Am,n - а2

2 m,n m,n

m,n Az m,n

bz m,n — q z m,n 1

Ax2 m,n Az

(Л.,п + 2а»,п )

2 2\ У m,n r~m,nj + Am,n

Az2 Ax2

L-±L m,n i-U m,n J

pm,n -а

C

—k

xm+0.5,n m,n * 2

Ax

2 (4 (m+1,n Лm-1,n)+ 2 (

А m+1,n Am-1,n

C ——k

x m-0.5,n m,n 2

Ax

)+ 2А"+0.5,п + Л.+0.5,п ),

2 (4 (m+1,n Лm-1,n)+ 2 (Am+1,n Am-1,n) 2а.-0.5,п Л.-0.5,п ),

C „<- — k

xm,n+0.5 m,n

1 , — k 1 ,

2 Am,n+0.5 , Cxm,n-0.5 m,n » 2 Am,n-0.5 -

Az m,n Az m,n

— km

1

Cym+0.5,n km,n » 2 Am+0.5,n , Cym-0.5,n km,n » 2 Am-0.5,n '

Ax m,n Ax m,n

— km

1

—k

^ym,n+0.5 lvm,n » 2 Am,n+0.5 , Cym,n-0.5 km,n » 2 Am,n-0.5 r

Az m,n Az m,n

=k

=k

^zm+0.5,n lvm,n l 2 Am+0.5,n , Czm-0.5,n km,n » 2 Am-0.5,n '

Ax m , n Ax m , n

=k

d — 1 k

xm+1,n 4 m,n

d = 1 k

xm,n+1 4 m,n

d = 1 k

z m+1,n 4 m,n

d = 1 k

z m,n+1 4 m,n

km,n a _2 (4 C^m,n+1 Лm,n-1 )+2 (Am,n+1 Am,n-1 )+ 2Am,n+0.5 + Л.,п+0.5 )

m,n

2 (4 (m,n+1 Лm,n-1)+ 2 (Am,n+1 Am,n-1) 2Am,n-0.5 Л.,п-0.5 ),

А +1 - А 1 ), d 1 — -d

,T~m,n+1 r~m,n-1/r’ xm—1,n x

z m,n+0.5 m,n 2

Az m,n

C ——k

z m,n-0.5 m,n 2

Az

fm

fm

— 1 k

m+0.5,n+0.5 4 m,n

Ax m. Az ,n m,n

1

Ax m. Az ,n m,n

1

Ax m. Az ,n m,n

1

Ax m. Az n m,n

1

1--------(л +1 — Л 1 ), d 1 —-d

Az m+1,n m-1,n x m,n-1

m,n

1--------(л 1 — Л 1), d 1 —-d

Az m,n+1 m,n-1 z m-1,n

m,n

A_ (Am+1,n Am-1,n ) d.

m,n

A (m+0.5,n+0.5 + Am+0.5,n+0.5 )

x m,n+1 *

z m,n-1 dz m,n+1 -

Axm,nAzm,n

m,n m,n

— 1 k

m-0.5,n-0.5 4 m,n

Ax Az

m,n m,n

m+0.5, n+0.5 + Am+0.5,n+0.5 m—0.5, n—0.5 + Am-0.5,n-0.5

),

m,n

m,n

1

1

1

1

/т

—-1к

т+0.5,п-0.5 4 т,п

■V-(Л

Ахт,пАгт,п

т,п т,п

/т

—-1к

т-0.5,п+0.5 4 т,п

+

■V-(Л

Ахт,пАгт,п

т,п т,п

'т+0.5,п-0.5 + Цт+0.5,п-0.5 , 'т-0.5,п+0.5 + Цт-0.5,п+0.5 ,

7 (Мп

т+0.5,п >*'*'т+0.5,п 1 4 УЛ*т+1,п Цт—1,п}^ ^1т-0.5,п ^*т-0.5,п 4 УЛ*т+1,п ^т-1,п/

1 _ ,, ^ _ 1

-1,п ) , Цт

Цт+0.5,п Ц

М т, п+0.5 Цт,п+0.5 1 4 т,п+1 т,п-1 т,п-0.5 т,п-0.5 4 т,п+1 Л* т,п

+ Т (Мт,п+1 - Мт,п-1 ), М т,

— М„

— М„

4 (цт+1,п Мт-1,п ) ,

(цт,п+1 Мт,п-1 ) ,

Р т,п

м2

р»

кт,п Рт,пЧ .

Использование данных соотношений позволяет рассчитывать особенности процесса распространения сейсмических колебаний с высокой степенью адекватности. Для примера на рис. 1 приведены для сравнения теоретические и реальные сейсмограммы, полученные в ходе выполнения исследовательских работ на участке шахтного поля ш. Красноармейская-Западная.

Для двухмерной задачи времена расчетов сравнительно малы, и шаг решетки можно положить независящим от направления и одинаковым по всей модели (ДА = тт{Дхтп, Дут п}). Тогда часть коэффициентов упростится:

( ( + 3М ) >

т,п т,п 2

Ьхт,п — Ьг т,п — Ч 2 - 2 . , 2 Рт,п - а

АА

т,п

V у

1

' у т,п ^1 ‘ к 1 2 ^ т,п^~т,п

V аа у

схт+0.5,п к т,п (4 (Лт+1,п Лт-1,п )+ 2 (,Мт+1,п Цт-1,п )+ 2Цт+0.5,п + Лт+0.5,п ^

сх т-0.5,п к т,п (4 (Лт+1,п Лт-1,п )+2 (цт+1,п Цт-1,п ) 2Цт-0.5,п Лт-0.5,п ) ,

с — с — к /М с — с — к /М

хт,п+0.5 у т,п+0.5 т,пН*т,п+0.5 > хт,п-0.5 ут,п-0.5 т,пН'т,п-0.5 >

с — с — к /М с — с — к /М

у т+0.5,п г т+0.5,п т,п т+0.5,п у т-0.5,п г т-0.5,п т,п т-0.5,п

с^ т,п+0.5 кт,п (4 (Лт,п+1 Лт,п-1 )+ 2 (,Мт,п+1 Цт,п-1 )+ 2Цт,п+0.5 + Лт,п+0.5 ) ,

сг т,п-0.5 к т,п (4 (Лт,п+1 Лт,п-1 )+ 2 (цт,п+1 Цт,п-1 ) 2Цт,п-0.5 Лт,п-0.5 ),

^хт+1,п 4 к т,п (цт,п+1 Цт,п-1 ) , ^хт-1,п ^х т+1,п ,

^хт,п+1 4 к т,п (Лт+1,п Лт-1,п ), ^хт,п—\ ^х m,n+1,

т+1,п 4 к т,п (Лт,п+1 Лт,п-1 ) , т-1,п т+1,п ,

т,п+1 4 к т,п (цт+1,п Цт—1,п ) , т,п-1 т,п+1 ,

/х т+0.5,п+0.5 /г т+0.5,п+0.5 4 к т,п (Лт+0.5,п+0.5 + Цт+0.5,п+0.5 ) ,

/х т-0.5,п-0.5 /г т-0.5,п-0.5 4 к т,п (Лт-0.5,п-0.5 + Цт-0.5,п-0.5 ) ,

100

200

І, МС

у—-Д^ЛЛД

---Vі-Ал

12

+М? . ..

Ад ^-4д

у \АА/-

Алл '

АДДлДл*- і

^Л/лЛ^лЛ 'л/*^А/ УЛ/Л**-*-

Сейсмограммы, полученные в ходе выполнения исследовательских

работ на участке шахтного поля ш. Красноармейская-Западная.

Вверху - реальные, внизу - полученные с помощью математического моделирования.

I*

= Л

хт+0.5,«-0.5 ./ гт+0.5,«-0.5

7кт,п(л,

т,п\ т+0.5,п-0.5 г^т+0.5

I*

= Л

= = - і к

х т-0.5,п+0.5 */ г т-0.5,п+0.5 4 т,п

Гл +'

т-0.5,п+0.5 “ ^т-0.5,п+0.5 У

12

_ _ і

V , 1 Л|^Луг'

у» -

- V у

- ■ "А *

■ - — у V/ 4

, Д Л Ж І

. ^ Дь V

- ^ і

V V у - —* А

- -"'Л

Остальные коэффициенты будут иметь тот же вид.

Необходимо отметить тот факт, что эти соотношения имеют самый общий вид. Какую бы физико-

математическую модель среды мы не использовали, слагаемые соотношений в рамках трехслойного конечноразностного представления будут те же. Изменять вид (или обнуляться) могут только коэффициенты при смещениях (например, для однородной среды все й^тп равны 0). В рамках данного анализа важно, что коэффициенты представляют собой линейные комбинации ^ и X в центральном узле и на половине расстояния до соседних. Поэтому изменение значений данных коэффициентов при учете дополнительных характеристик углевмещающей толщи и ее состояния можно представить в виде условной замены ^ и X на ^эфф и Хэфф:

к =

т,п

АИ

-Рт,пЯ .

100

200

І, МС

1

1

Иэфф = И + Ди,

Кфф = Л + Дх.

где Ди и ДЛ рассматриваются как дополнительные слагаемые, вносимые за счет учета состояния горного массива и характеристик пород.

Такого подход обоснован и с физической точки зрения [1, 4, 7, 8]. Например, в монографии [8] на основе анализа результатов экспериментальных исследований, проведенных на ряде шахт Украины и России, представлены зависимости модуля сдвига и других параметров среды от коэффициента пористости, от величины горного давления и других парамет-

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

Одним из наиболее сложных вопросов является учет затухания колебаний. Оно обусловлено целым рядом факторов [1], среди которых важнейшими являются расхождение фронта волны, рассеяние волны на неоднородностях, поглощение упругих волн и др. Затухание колебаний за счет расхождения фронта волны при использовании МКР учитывается автоматически. Рассеяние волны на неоднородностях может быть промоделировано в случае, если размер неоднородности сравним с шагом дискретизации модели. В то же время моделирование поглощения волн представляет собой чрезвычайно сложную задачу.

Колебательный процесс выражается в периодических сжатиях, растяжениях частиц среды и смещениях их относительно друг друга. Определенная доля энергии колебаний при этом расходуется. Характер и параметры этого процесса существенно зависят как от типа волны, так и от физико-механических свойств среды. Он описывается коэффициентом поглощения в и логарифмическим декрементом затухания в, которые связаны между собой соотношением в = рь, где Ь - длина волны. При расчете компонент смещений волнового поля в заданной точке неоднородной среды учесть поглощение аналитически в общем случае практически невозможно в результате того, что колебания достигают источника по-разному, проходя различное расстояние по разным породам, претерпевая отражения, преломления и трансформации на границах раздела сред. В итоге, в настоящее время доминирует упрощенный подход, рассматривающий величину р как интегральный показатель, зависящий от параметров среды, от размеров нарушенной области, от соотношения размеров нарушенной и ненарушенной зон. Учет интегрального поглощения имеет «оценочный» характер и производится с помощью фильтрации теоретических сейсмограмм.

В рамках МКР можно предложить подход, который в явном виде учитывает реальные пути распространения колебаний и трансформации на границах раздела сред и на неоднородностях. Рассмотрим физический смысл соотношений (3)-(5). Можно считать, что они описывают математическую модель, в рамках которой любое смещение в центре элемента етп в момент времени р+1 представляется в виде суперпозиции различных типов слагаемых. Первый из них описывает волны, пришедшие за единичный временной интервал с восьми направлений (от центров соседних элементов). Вне зависимости от того, как долго эти волны реально были в пути, они прошли известный путь (расстояние между центрами элементов) по известным породам. Этот факт дает право в коэффициенты при смещениях в соседних с етп элементах ввести поправки на поглощение по правилу, показанному ниже на примере коэффициентов сх т+0.5,п и Лх т+1,п:

Сх т + 0.5,п = Сх т+ 0.5,п Х еХр(— "2 Рх т,п^Хт,п — ~ Рх т+1,п^Хт+1,п ) ,

<т+1,п = Лхт+1,п Х ехр(- ^ р^ Д^п - Т в т+1,п АХт+1,п ) ,

где рх т п и ръ тп представляют собой коэффициенты поглощения РУ и 8У волн соответственно в элементе ет,п с учетом пористости, трещиноватости, величины горного давления и других факторов. Во втором примере мы имеем дело с трансформацией волны, поэтому считаем, что по элементу ет+1п волна распространялась как БУ, а по ет+1п - как РУ. Записывать остальные коэффициенты не будем, дабы не загромождать рассмотрение. Второй тип слагаемых можно трактовать как остаточную часть от колебаний в предыдущие моменты времени в том же элементе (при расчете которых ранее уже учитывался фак-

тор поглощения). Зависимость в от пористости, трещиноватости и величины горного давления подробно рассмотрена в ряде работ [І, 5, В] и поэтому в рамках данной статьи не рассматривается.

Конечно, данный подход, как и общепринятый, имеет свои недостатки. Основную сложность вносит тот факт, что коэффициенты поглощения зависят от частоты. Эта зависимость достаточно хорошо изучена [І, 5, В] и, вообще говоря, имеет сложный нелинейный характер. Этот факт вносит ограничение на диапазон применимости данного подхода. Мы можем считать его корректным в случае, если в рамках конкретной решаемой задачи в пределах частотного диапазона исследуемой части теоретического сигнала данным фактором можно пренебречь. Например, опираясь на приведенные в работе [В] графики модельных зависимостей данного параметра от частоты, можно сделать вывод о том, что в пределах низких частот (до 200-250 Гц) они имеют практически линейный характер с относительно слабым возрастанием. Чем ниже частота исследуемого волнового пакета и уже его спектр, тем обоснованнее применимость данного метода в «чистом виде».

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

В настоящей работе рассмотрены конечно-разностные соотношения, позволяющие моделировать процесс распространения сейсмических колебаний в угленосной толще. Предложены способы учета напряженного состояния горного массива, пористости, трещиноватости угля и пород при математическом моделировании процесса распространения сейсмических колебаний в угленосной толще и учета поглощения сейсмических колебаний в среде при использовании МКР. Данный подход реализован в комплекте программ, который в настоящее время используется в институте УкрНИМИ для математического моделирования процессов распространения сейсмических колебаний в угленосной толще. Опыт применения МКР показал, что данный подход особенно актуален в сложных горно-геологических условиях и одинаково успешно может быть применен как при решении задач наземной, так и подземной сейсморазведки.

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

------------------------------------------------------------------ СПИСОК ЛИТЕРАТУРЫ

1. Азаров Н. Я., Яковлев Д. В. Сейсмоакустический метод прогноза горно-геологических условий эксплуатации угольных месторождений. - М.: Недра,І9ВВ. -І99 с.

2. Korn M. Stocl H. Reflection and Transmission of Love channel Waves at Coal Seam Discontinuitis Computed with A Finite-Difference Method. -J. Gejphis.,1982,50. p. І7І-І76.

3. Глухов А.А., Захаров В.Н., Рубан А.Д. Моделирование волнового поля в задачах шахтной сейс-

моразведки методом конечных разностей/Горный вестник, М., изд. ИГД Скочинского, І994, С.І6-ІВ

4. Анциферов А.В. Моделирование волнового поля в задачах шахтной сейсморазведки методом конечных разностей/ Збірник наукових праць №5 «Проблеми гірського тиску» 200І. С.5-І5.

5. Анциферов А.В. Теория и практика шахтной сейсморазведки. - Донецк: изд. «Алан», 2002, -3І2 с.

6. БреховскихЛ.М. Распространение волн в слоистых средах. - М.: Наука, І973.

7. Рубан А.Д., Захаров В.Н. Исследова-ние зон повышенного горного давления (ПГД) в

углепородных массивах сейсмоакустическим методом // Механика горных пород: Науч. Сообщ./ ННЦ ГП-ИГД им. А.А. Скочинского, - М. - І999. - Вып. З ІЗ. - С. 39-4В.

В. Захаров В.Н. Сейсмоакустическое прогнозирование и контроль состояния и свойств горных

пород при разработке угольных месторождений. - М.: ФГУП ННЦ ГП - ИГД им. А.А. Скочинского, 2002. -І72 с.

Коротко об авторах

Глухов А.А. - кандидат технических наук, заведующий отделом, Анциферов А.В. -доктор технических наук, директор, УкрНИМИ, Донецк, Украина.