Кормильцев В.В., Ратушняк А.Н. Моделирование геофизических полей при помощи объемных векторных интегральных уравнений

Подождите немного. Документ загружается.

ления, что позволяет представить аномальное магнитное поле в фор-

ме (13).

Однако в случае полупространства магнитное поле токов, расте-

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

следует вернуться к закону Био–Савара. Однако, как показал Стефа-

неску [6, c.315–317] , для точечного источника постоянного тока, ин-

тегрирование электрических токов по всему объему нижнего полупро-

странства можно заменить изменением конструкции сторонних токов,

добавив к ним полубесконечный кабель, направленный перпендику-

лярно границе раздела в сторону от нее. При определении магнитного

поля в воздухе кабель начинается в месте заземления, а при опреде-

лении в земле – в месте изображения заземления в границе раздела.

Если ось OZ направлена вверх и заземление является источником то-

ка, а не стоком, то ток в каждом из кабелей следует направить вниз.

Этот фиктивный ток является как бы продолжением тока в диполе.

Подробно метод изображений при вычислении магнитного поля токов

рассмотрен в [7,c.74–77]. Применительно к векторному потенциалу это

означает появление дополнительного слагаемого в составляющей

Поскольку каждый дипольный источник представляет собой два сбли-

женных точечных заземления, соображения Стефанеску можно при-

менить к нашему случаю.

Обозначим момент диполя

M, а момент отраженного диполя M

Они различаются тем, что

=–M

. Добавка в П

составляет при z≤0

πΠ

ζζ

rz rrz

−+

⋅

−+()[(

)]

что вместе с прежним значением векторного потенциала диполя дает

новое значение [8]

111

⋅

−+r

rr z

()

[( )]

причем

=− =

⋅

⎡

⎣

⎢

⎤

⎦

⎥

div П Mr M r

σπσ

что доказывает правильность выкладок. Двойной знак в формулах по-

зволяет учесть свойства верхнего полупространства. Применительно

к электрическому и магнитному полю тока верхний знак соответствует

изолятору (воздуху), нижний – идеальному проводнику. Выражение

для

проще всего получить путем предельного перехода от наклон-

ного заземленного прямолинейного кабеля к диполю. С учетом этих

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

при

z≤0, приравняв M к i

пa

, получим

HrotП rot

ik ir

rr z

dV==

+⋅

−+

⎧

⎨

⎩

⎫

⎬

⎭

∫

πζ

п a п az п a1

()

[( )]

(14)

Конкретный вид

па

и определяющих его сторонних токов i

зави-

сит от сущности процесса. Далее будут рассмотрены два явления в

горной породе, содержащей в своих порах раствор электролита, элек-

тронейтральность которого нарушена, а подвижности ионов изменены

за счет адсорбции ионов одного знака минеральным скелетом. Это

сторонние токи течения Дарси и диффузии. Будет рассмотрен и са-

мый общий случай, когда в отношении стороннего тока делают един-

ственное предположение, а именно, что он вызван предшествующим

протеканием электрического тока и пропорционален его плотности

=–

i. Коэффициент пропорциональности называют поляризуемо-

стью, а само явление – вызванной поляризацией.

Однако прежде чем обратиться к электрическим явлениям, рас-

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

уравнениями.

1.3. Аналогичные интегральные уравнения.

Имеется шесть физических задач, приводящих к уравнениям эл-

липтического типа. Кроме уже рассмотренной задачи о растекании по-

стоянного электрического тока к уравнениям Лапласа и Пуассона при-

водят задачи электростатики, магнитостатики, стационарные задачи

теплопроводности, диффузии и течения Дарси. Их решения аналогич-

ны, аналогичны и условия сопряжения на границах разрыва физиче-

ских свойств. Следующая схема, позаимствованная из [3], иллюстри-

рует подобие перечисленных полей:

Электрический по-

тенциал

Удельная электро-

проводность

Плотность электриче-

ского тока

E=–

gradU

Электростатический

потенциал

Диэлектрическая

проницаемость

Электрическая индукция

E=–

gradU

Магнитный потен-

циал

Магнитная прони-

цаемость

Магнитная индукция

H=–

gradU

Давление P

Отношение гидрав-

лической прони-

цаемости к динами-

ческой вязкости

Скорость течения Дарси

v=–c/

gradP

флюида

Температура T

Коэффициент теп-

лопроводности

Плотность теплового

потока

q=–

gradT

Концентрация C

Коэффициент диф-

фузии

Плотность потока веще-

ства

q=–D gradC

По аналогии между напряженностью электрического поля токов и

магнитного поля магнетиков запишем для напряженности магнитного

поля

H неоднородного магнетика в немагнитной среде (

=1, I

=0)

[5]

[]

HH grad HI

=− −+

∫

()

(15)

где

– напряженность внешнего магнитного поля, например земного;

– магнитная проницаемость неоднородного магнетика; I

– его неод-

нородная остаточная намагниченность. Для слабомагнитных объектов

и в слабых полях

μ−

=κ

, где

– магнитная восприимчивость; для

ферромагнетиков в сильных полях значения

можно уточнять на каж-

дом шаге итерации, пользуясь вычисленной напряженностью внут-

реннего поля и кривой намагничивания.

Для течения Дарси [5,8]

grad grad grad

Kr K r

=−

⋅

⎡

⎣

⎢

⎤

⎦

⎥

∫

Kij=−

⎛

⎝

⎜

⎞

⎠

⎟

⎛

⎝

⎜

⎞

⎠

⎟

PPP

∂

∂ξ

∂

∂η

∂

∂ζ

Kij

1=−

⎛

⎝

⎜

⎞

⎠

⎟

+−

⎛

⎝

⎜

⎞

⎠

⎟

PPP

∂

∂ξ

∂

∂η

∂

∂ζ

(16)

где

– давление, развиваемое источниками флюида в однородном

по гидравлической проницаемости пористом полупроcтранстве, c и

– проницаемости неоднородности и полупространcтва. Верхний знак в

уравнении для случая, когда нижнее полупространство гидравлически

изолировано от верхнего и

∂P/∂z=0 при z=0, нижний – для случая, ко-

гда дневная поверхность является поверхностью высачивания и

P=0

при

z=0.

По аналогии между плотностью постоянного электрического тока

и плотностью стационарного потока тепла запишем интегральное

уравнение для градиента температуры в нижнем полупространстве,

если на его поверхности поддерживается постоянная температура, в

виде [5,9]

grad grad grad

Gr G r

=−

⋅

−

⋅

⎡

⎣

⎢

⎤

⎦

⎥

∫

Gij=−

⎛

⎝

⎜

⎞

⎠

⎟

⎛

⎝

⎜

⎞

⎠

⎟

∂

∂ξ

∂

∂η

∂

∂ζ

TTT

Gij

1=−

⎛

⎝

⎜

⎞

⎠

⎟

+−

⎛

⎝

⎜

⎞

⎠

⎟

∂

∂ξ

∂

∂η

∂

∂ζ

TTT

(17)

где

T – температура; q=–

gradT

=–

∂T

/∂z – плотность теплового

потока, идущего снизу к дневной поверхности;

– коэффициенты

теплопроводности полупространства и неоднородности.

Рассмотрим искажения стационарного потока диффун-

дирующего растворенного вещества или электролита, вызванные не-

однородностью пористой среды. Плотность потока

q удовлетворяет

закону Фика

q=–mDgradC , где m – коэффициент пористости, D – ко-

эффициент диффузии,

C – концентрация растворенного вещества,

например, газа или электролита. В случае растворенного электролита

на движение ионов влияет не только осмотическое давление, но и

электрическое поле. В пористой среде электронейтральность раство-

ра нарушена, что изменяет подвижности аниона и катиона, а также ко-

эффициент диффузии электролита [10]. Коэффициент диффузии

D в

порах зависит от формы и размера пор, не связан линейно с коэффи-

циентом пористости

m и в общем случае является функцией коорди-

нат, независящей от распределения пористости. В связи с этим необ-

ходимы пояснения относительно уравнения неразрывности и условий

сопряжения на границах сред. Поток вещества неразрывен, но кон-

центрация электролита в порах испытывает скачок на границе двух

пористых сред, поддерживаемый скачком мембранного потенциала.

Поэтому

div(mDgradC )=0, C

≠C

, m

∂C

/∂n= m

∂C

/∂n.

Однако отношения

C/C

и D/D

, где C

– концентрация электроли-

та в свободном растворе, равновесном с пористой средой, и

– ко-

эффициент диффузии в растворителе, мало отличаются от единицы

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

рим с толщиной двойного слоя [11]. Поэтому будем считать коэффи-

циент диффузии всюду постоянным, а концентрацию непрерывной.

Тогда

div( )mC grad = 0

∂

= ,

(18)

∂

21 2

1−=−

⎛

⎝

⎜

⎞

⎠

⎟

и по аналогии с (8) приходим к интегральному уравнению для гради-

ента концентрации в безграничной пористой среде с неоднородным

включением:

grad grad grad grad

=− −

⎛

⎝

⎜

⎞

⎠

⎟

∫

(19)

где

gradC

– градиент концентрации диффундирующего вещества в

отсутствие неоднородности.

1.4. Родственные интегральные уравнения

Выпишем интегральные уравнения, родственные уравнению (8)

для электрического поля, порождаемого объемно распределенными

сторонними токами. Рассмотрим проводящее нижнее полупространст-

во, в котором после пропускания постоянного тока установилась ста-

ционарная вызванная поляризация. Положим в (2), что

cn111

ση

∂

cn222

ση

∂

, где

– установившиеся значения поляризуе-

мости полупространства и неоднородности. Тогда вместо (3) имеем

∂

ση

∂

21 22

−=−

−

⎛

⎝

⎜

⎞

⎠

⎟

()

Повторяя рассуждения и выкладки, приводящие к выражению

(11), получаем для суммарного поля, опустив индекс 2 для неодно-

родности

grad

Er E r

−

⎛

⎝

⎜

⎞

⎠

⎟

−

⎛

⎝

⎜

⎞

⎠

⎟

⋅

⎛

⎝

⎜

⎞

⎠

⎟

∫

πσ η π

ση

()

(20)

Определяя напряженность поля вызванной поляризации, необ-

ходимо решить это интегральное уравнение еще раз, положив

=0, и найти разность.

Пусть в проводящем нижнем полупространстве в результате те-

чения раствора электролита в порах горной породы возникают потен-

циалы течения (фильтрационные) и электрическая напряженность те-

чения Дарси

Согласно [6, c.239–240], положим

LgradP , где P – давление,

L – коэффициент потенциала течения. Тогда вместо (3) получим, что

∂

21 2

−=−

⎛

⎝

⎜

⎞

⎠

⎟

+−

=−

⎛

⎝

⎜

⎞

⎠

⎟

+−

⎛

⎝

⎜

⎞

⎠

⎟

(21)

Последнее равенство получается после учета условия непре-

рывности нормальной к поверхности

S составляющей скорости тече-

ния Дарси

11 2 2

2nn

=− =− =

∂

∂μ

∂

где

и c

– проницаемости неоднородного объема и среды,

– дина-

мическая вязкость флюида. Используя (21), получаем на основании

выкладок, приводящих к уравнениям (8), (11), что

E E grad

Wr W r

=−

⋅

⎡

⎣

⎢

⎤

⎦

⎥

∫

Egra

−

⋅

()

Wijk ij=−

⎛

⎝

⎜

⎞

⎠

⎟

++ + −

⎛

⎝

⎜

⎞

⎠

⎟

⎛

⎝

⎜

⎞

⎠

⎟

∂

∂ξ

∂

∂η

∂

∂ζ

ξηζ

1 EEE LL

PPP

()

Wijk ij

1=−

⎛

⎝

⎜

⎞

⎠

⎟

+− + −

⎛

⎝

⎜

⎞

⎠

⎟

+−

⎛

⎝

⎜

⎞

⎠

⎟

∂

∂ξ

∂

∂η

∂

∂ζ

ξηζ

EEE LL

PPP

(22)

Верхний знак “плюс” в выражении (22) соответствует случаю, когда

нижнее полупространство изолировано от верхнего тонким слоем не-

проницаемых пород (

∂P/ ∂z =0 и E

=0 при z =0). Знак минус соответст-

вует случаю, когда дневная поверхность является поверхностью вы-

сачивания (

P=0). При этом она одновременно является поверхностью

равного потенциала (

=0 при z=0).

Для расчетов напряженности электрического поля течения Дар-

си, согласно (22), необходимо предварительно рассчитать фильтра-

ционную задачу (16) при совпадающих внешних границах неоднород-

ного объема

V. В этом случае матрица внутренних значений gradP ,

полученная в результате решения уравнения (16), без каких-либо из-

менений целиком используется в уравнении (22).

Рассмотрим напряженность электрического поля диффузии

электролита, положив

gradC, где

=RT(v

–v

), R – универсальная

газовая постоянная,

T – абсолютная температура, v

и v

– скорости

движения аниона и катиона в единичном электрическом поле [12]. Из-

менения

в пространстве вызваны изменением величины просвета

между заряженными стенками пор. В данном случае изменениями

нельзя пренебречь, как пренебрегли ранее изменениями коэффици-

ента диффузии, поскольку это существенная сторона рассматривае-

мого явления. Подставляя значение

в (3) и учитывая (18), имеем

∂

∂σ

ββ

∂

21 2

222

211

−=−

⎛

⎝

⎜

⎞

⎠

⎟

+−

=−

⎛

⎝

⎜

⎞

⎠

⎟

+−().

Тогда интегральное уравнение для электрического поля диффузии

электролита в пористой среде может быть записано в виде

E E grad E grad

=− −

⎛

⎝

⎜

⎞

⎠

⎟

+−⋅

⎡

⎣

⎢

⎤

⎦

⎥

∫

σσ

ββ

()

(23)

где

Egra

=−

– напряженность первичного поля диффузии в

отсутствие неоднородности,

gradC

– внутренние значения градиента

концентрации, полученные при решении уравнения (19).

Выпишем интегральное уравнение для электроосмоса, родст-

венное уравнению (16). При наложении электрического поля

E в по-

ристой среде возникает электроосмотическое течение, скорость и гра-

диент давления которого равны

v=–Li=–

LE, gradP=

E/ c , где c –

проницаемость,

– динамическая вязкость флюида. Согласно соот-

ношениям Онсагера, коэффициент, связывающий плотность электри-

ческого тока со скоростью электроосмотического течения, совпадает с

коэффициентом потенциала течения L [11]. Используя аналогию меж-

ду течением несжимаемой жидкости в пористой среде и электриче-

ским током, между давлением и потенциалом, вместо (2) и (3) записы-

ваем

−− =− −

11 1

22 2

∂

μσ

LL E

n21 2

2122

1−=−

⎛

⎝

⎜

⎞

⎠

⎟

−()

с учетом того, что

. Тогда для градиента давления при

электроосмосе

grad grad grad

Gr G r

=−

⋅

⎡

⎣

⎢

⎤

⎦

⎥

∫

Gijk=−

⎛

⎝

⎜

⎞

⎠

⎟

⎛

⎝

⎜

⎞

⎠

⎟

+−

PPP

∂

∂ξ

∂

∂η

∂

∂ζ

μσ

()

Gijk

1=−

⎛

⎝

⎜

⎞

⎠

⎟

+−

⎛

⎝

⎜

⎞

⎠

⎟

+−

PPP

∂

∂ξ

∂

∂η

∂

∂ζ

μσ

()

(24)

grad EP

Верхний знак в уравнении (24) – для случая, когда нижнее полупро-

странство гидравлически изолировано от верхнего и ∂P/∂z =0 при z=0,

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

стью высачивания и P =0 при z =0. Там, где разность

gradP– gradP

≠0,

возникают локальные области напряжения или разгрузки. Наиболее

вероятная причина интенсивного электроосмоса – электротеллуриче-

ская напряженность

i+E

j, cвязанная с магнитными пульсация-

ми во время общей геомагнитной бури. Если объем V отличается по

электропроводности от вмещающих пород, предварительно необхо-

димо решить уравнение

E E grad E

=− −

⎛

⎝

⎜

⎞

⎠

⎟

⎛

⎝

⎜

⎞

⎠

⎟

∫

(25)

и использовать полученную матрицу значений внутреннего поля E

, E

в уравнении (24).

1.5. Интегральные уравнения для градиентов температуры и кон-

центрации в фильтрующей среде.

Рассмотрим векторное интегральное уравнение для градиента

температуры в случае комбинированного переноса тепла кондуктив-

ным и конвективным способами [13]. Конвективный перенос осущест-

вляется за счет фильтрации флюида в проницаемой влагонасыщен-

ной среде по закону Дарси со скоростью

v=–c/

gradP, причем эта ско-

рость невелика, так что температура флюида и минерального скелета

в каждой точке успевает стать одинаковой. Среда представляет собой

однородное пространство или полупространство с постоянными ко-

эффициентами гидравлической проницаемости

и теплопроводности

Для случая стационарной температуры при столь медленной

фильтрации, что температура T твердой и жидкой фаз одинакова, из

закона сохранения энергии [14] имеем для дивергенции полного пото-

ка тепла, что

()

[]

div

λρ

⋅− +grad vTcTH

(26)

Здесь

– коэффициент теплопроводности, Вт/м/

;

– плотность

флюида, кг/м

; c

– удельная массовая теплоемкость флюида,

Дж/кг/ °К; H – тепловыделение в единице объема, Вт/м

, которое мо-

жет иметь радиогенное, химическое или биохимическое происхожде-

ние;

v=–c/

gradP – скорость течения Дарси, м/с; c – гидравлическая

проницаемость, м

;

– динамическая вязкость флюида, Па⋅c; P – дав-

ление в порах, Па. Отсюда в случае несжимаемой жидкости

T+ −

grad

=−

λχλ

(27)

поскольку div(

v )=0. Здесь

⋅

– коэффициент, подобный темпе-

ратуропроводности, м

/с.

Cначала рассмотрим безграничную среду, в которой находится

одно включение объема V, имея в виду, что впоследствии интеграль-

ное уравнение может быть обобщено на случай нескольких объемов

непосредственно. Напомним, что в среде

grad

=0 и во включении

grad

≠0. Представим температуру в виде T=T

, где T

– темпера-

тура в среде без включения, удовлетворяющая уравнению

=−

grad

χλ

(28)

Это уравнение подобно уравнению Пуассона, где в правой части стоит

объемная концентрация источников тепла. Решением его является

интеграл Пуассона

HdV

=− −

⎛

⎝

⎜

⎞

⎠

⎟

∫

πχ λ

grad

=(x–

)

+(y–

)

+(z–

)

, dV=d

ξ⋅

η⋅

Интегрирование ведется по всему безграничному пространству. Инте-

грал Пуассона является интегро-дифференциальным уравнением от-

носительно температуры. Дифференцируя его по координатам точки

наблюдения, получаем интегральное уравнение относительно гради-

ента температуры. Решение уравнения (28) может быть получено так-

же методом разделения переменных или еще какими-либо способами.

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

температуру T

полагают всюду заданной.

Приведем здесь несколько иную процедуру вывода выражения

для аномальной температуры, нежели в разделе 1.1, имея ввиду, что

конечные формулы тождественны. Использовав аналогию между тем-

пературой и электрическим потенциалом, запишем, что [4]

=+ +

∫∫

λπ

(29)

где

=–divF – избыточная концентрация объемных источников тепла

во включении, Вт/м

;

– концентрация поверхностных источников

тепла, Вт/м

; S – поверхность включения, м

; F, F

– вектор поляриза-

ции и его нормальная составляющая, Вт/м

Второе и третье слагаемые соответствуют аномальной тем-

пературе T

, вызванной включением. Избыточная концентрация объ-

емных источников тепла равна разности правых частей уравнений (27)

и (28):

λχχλ

=−+−

grad

gradTT T

(30)

Плотность поверхностных источников пропорциональна скачку

нормальной составляющей градиента температуры между внутренней

и внешней сторонами граничной поверхности

(

)

(

)

∂

=−

+−

Здесь под T

(+)

и T

(–)

можно понимать как суммарные T, так и аномаль-

ные T

температуры. Граничными условиями являются непрерывность

температуры T

(+)

(–)

и нормальной составляющей полного потока те-

пла. Согласно (26),

(

)

(

)

(

)

(

)

(

)

∂

ρλ

∂

ffn ffn

−

−=−

()

Очевидно, что движение флюида не создает поверхностных источни-

ков тепла, поскольку фильтрационную задачу решают при граничном

условии v

(+)

(–)

. Поэтому

(

)

(

)

(

)

∂

+−

−=−

⎛

⎝

⎜

⎞

⎠

⎟

Это выражение является нормальной составляющей F

вектора поля-

ризации

F единицы объема включения

Fgr=−

⎛

⎝

⎜

⎞

⎠

⎟

Tad

(31)

Подставляя (30) и (31) в (29) и заменяя в полученном выражении по-

верхностный интеграл на объемный по формуле Гаусса

FdS

div

⎛

⎝

⎜

⎞

⎠

⎟

∫∫

имеем после выполнения операции дивергенции над функцией

F/r,

замены

T на правую часть выражения (27) и приведения подобных

членов, что

TT T

TTHH

=+ −

⎛

⎝

⎜

⎞

⎠

⎟

⎛

⎝

⎜

⎞

⎠

⎟

−

⎛

⎝

⎜

⎞

⎠

⎟

∫

πχ λ

grad grad

vgrad v grad

(32)

Здесь индекс M означает дифференцирование по координатам

Bторое слагаемое представляет собой дипольно поляризованную, а

третье – источниковую часть поля аномальных температур. Совер-