CC BY
71
Геофизика
УДК 550.831.01
МЕТОДЫ ТЕОРИИ ПОТЕНЦИАЛА ПРИ РЕШЕНИИ ПРЯМЫХ ЗАДАЧ ГРАВИМЕТРИИ И ГЕОДИНАМИКИ ТРЕХМЕРНЫХ НЕОДНОРОДНЫХ СРЕД
Ю.В. Пятаков, В.И. Исаев*, В.Ю. Косыгин**
Воронежский государственный университет инженерных технологий *Томский политехнический университет “Вычислительный центр ДВО РАН, г. Хабаровск E-mail: pyatakovjv@mail.ru
Приведены алгоритмы решения задачи определения составляющих гравитационного и геодинамического полей для трехмерных неоднородных сред. Для аппроксимации плотностной и реологической структуры среды используются типовые элементы -вертикальные треугольные призмы с произвольными верхним и нижним основаниями. Получено новое аналитическое решение прямой задачи гравиметрии для типового элемента с плотностью, меняющейся с глубиной по линейному закону. Выполнена математическая постановка и приведено общее решение задачи определения напряжений и мгновенных скоростей смещения неоднородной вязкой среды под действием гравитационного поля Земли. Решение определено с использованием гидродинамических потенциалов объемного, простого идвойного слоя. Показано, что для численного расчета таких потенциалов оптимально использовать теорию и типовую технику решения прямых задач гравитационного потенциала. Устойчивость, точность и быстродействие разработанных алгоритмов демонстрируется расчетами тестовых примеров.
Ключевые слова:
Теория потенциала, прямые задачи гравиметрии и геодинамики, алгоритмы.
Key words:
Potential theory, direct gravity and geodynamics problem, algorithms.
Введение
Фундаментальной основой методов гравиметрии является теория потенциала, первоначально понимаемая как учение о свойствах сил, действующих по закону всемирного тяготения. В дальнейшем в работах К.Ф. Гаусса, Ж.Л. Лагранжа, П.С. Лапласа, Д. Грина, Ж.А. Пуанкаре, А.М. Ляпунова методы теории потенциала были распространены не только на решение задач теории тяготения, но и на решение широкого круга задач математической физики.
В настоящее время теория потенциала является мощным и универсальным инструментом, широко используемым при решении задач интерпретации геофизических полей. Одним из примеров такого использования является известное соотношение Пуассона, устанавливающее связь между решениями прямых задач магниторазведки и гравиразведки. Это позволяет использовать богатый арсенал средств, накопленных в гравиметрии, для расчета и интерпретации магнитных аномалий.
В настоящей статье рассматриваются вопросы построения алгоритмов решения прямых задач гравиметрии и геодинамики. Под прямой задачей геодинамики здесь понимается задача расчета параметров напряженно-деформированного состояния среды, обусловленного неоднородностью ее плотностного строения.
В основу решения задач положены методы теории потенциала, что позволяет создавать эффективные вычислительные алгоритмы на базе типовой методики решения прямых задач гравиметрии, изложенной в работах [1-4].
Решение прямой задачи гравиметрии
для трёхмерных блоково-градиентно-слоистых сред
Как отмечено в работах [4-7], при плотностном моделировании геологических структур в гравитационном поле их необходимо сводить к трехмерным моделям, а плотностную параметризацию осуществлять с учётом изменения плотности как в латеральном направлении, так и в направлении от кровли к подошве.
В качестве расчетного элемента в этом случае удобно использовать вертикальную треугольную призму (рис. 1), плотность которой записывается в виде линейной функции отглубины p(^ъ)=po+k^ъ. Здесь и далее, для сокращения объёма записи, будем пользоваться взамен буквенной (используемой вработах [5-7]) числовой индексацией координат (используемой в работах [2, 3]), полагая х=хь y=X2, z=Xз - в обозначении координат точек расчета и 4=4, Ц=4, 0=4 - в обозначении координат переменных интегрирования.
Приводим решение прямой задачи гравиметрии (1), выполненное в соответствие с методикой работ В.Н. Страхова [2, 3]. Особенность этого решения состоит в том, что оно, в отличие от [2, 3], дано для тел с линейным законом изменения плотности.
Введем следующие обозначения: пусть ^и 52 -соответственно верхняя и нижняя грани призмы, а Sз, и $5 — её боковые грани. Тогда аналитическое выражение для вертикальной составляющей гравитационного потенциала можно записать в виде:
К, О) = -/Е[с'(х) +кг (х)]г3} + /'^ (х)- (2)
1=1
Здесь /кУ^) - гравитационный потенциал от призмы с постоянной плотностью к; V (х) определяется по формуле:
V(х) =Е Ф'; (х)/2, (3)
N і. .Ш(2) - Л(1))
V (х) = 2І іагіЬ ЛіЩ | аіап ' ^т т ^ , (4)
І=1
£(2) £(1) + і2 5
г(х) = £«„і = 1,2,...,5;
(5)
І=1
где
Рис. 1. К решению прямой задачи гравиметрии. Аппроксимирующая вертикальная треугольная призма. Условные обозначения и пояснения в тексте
В этом случае выражение для вертикальной составляющей гравитационного потенциала от призмы примет следующий вид:
V*, (х) = / Ш(Ро + к^)(4 - *3)^^^ (1)
В
где В - вертикальная треугольная призма,
р з)7;
/ - константа гравитации.
Об аналитических решениях прямой задачи гравиметрии для выделенного аппроксимирующего элемента сообщалось. В работах [7-10] приводятся принципиально идентичные, но отличающиеся по форме аналитические выражения для составляющих гравитационного потенциала от тел переменной плотности.
В.Н. Страховым в работе [2] был предложен новый подход к решению прямых задач гравиметрии для многогранных тел, который, как показали сопоставления, проведенные в работе [3] для тел с постоянной плотностью, позволяет получать вычислительные алгоритмы, превосходящие аналогичные решения по таким показателям, как точность, быстродействие и устойчивость.
А-, Д, €ь В1 - коэффициенты уравнения 1-й грани призмы Ах +Вх2+Схз +В=0, определяемые по значениям координат вершин грани; N - количество вершин ¡-й грани призмы; \Ь/ - длинау-й стороны грани (/=1,2,...,^); Яц - расстояние от расчетной точки до у-й вершины грани ^(Д,м+1=Д,1); х$, х/, х3/ - координаты у-й вершины грани ^.
Отметим, что формулы (2)-(5), также как и соответствующие формулы вработах [2, 3], являются едиными для расчета как внешнего (по отношению к телу В), так и внутреннего полей. Это позволяет использовать полученные формулы для обработки данных как наземных, так и скважинных измерений.
і=і
Таблица 1. Результаты вычисления вертикальной составляющей гравитационного потенциала в точках профиля, проходящего по ребру (профиль I) и внутри (профиль II) тела
Профиль Координаты точки расчета, км Расчет по формулам (2)-(5), 10-5м/с2 Расчет методом численного интегрирования (в Майлса^, 10-5м/с2 Расчет из работы [7], 10-5м/с2 Номер точки
Х Х2 Х3
4 0 1 1,06048024-102 - 1,06048024-102 1
8 0 2 6,78016491-10’ - 6.78016500-101 2
I 8,001 0 2,0025 6,77264793-10’ - 6,77264793-Ю1 3
12 0 3 2,15708388-Ю1 2,15708388-Ю1 2,15708397-Ю1 4
16 0 4 8,92005655-Ю1 8,92005655-Ю1 8,92005712-Ю1 5
Время расчета, с 3,9-10-5 1,8-10-1 7,3
-6 3 10 -1,35443390-10° -1,35443390-100 -1,35443390-100 1
1,999 3 10 5,77334110-10° - 5,77334110-100 2
II 2 3 10 5,77954460-100 - 5,77954460-100 3
2,001 3 10 5,78574882-100 - 5,78574882-100 4
4 3 1 1,46939534-Ю1 - 1,46939534-Ю1 5
Время расчета, с 4,5-10-5 9,0-10-2 4,5
Примечание. Здесь ив табл. 2 подчеркнуты расхождения с расчетом по аналитическим выражениям (2)-(5). Алгоритм Ма^сад, реализующий решение прямой задачи численным методом в точках, лежащих на ребре призмы, в точке 3 продолжения ребра, во внутренних точках тела и в точке 2 профиля II, «не работает».
Тестирование алгоритма решения
прямой задачи гравиметрии
Устойчивость, точность и быстродействие алгоритма решения прямой задачи гравиметрии демонстрируются расчетами на трех прямолинейных профилях, различно расположенных по отношению к телу. В качестве модельного примера используется пример из работы [7]. В этой работе тестовые расчеты аналитических решений проверялись в сопоставлении с расчетами по «прямому» численному алгоритму, реализующему метод Таль-вани [11].
Координатное описание тела, представленного на рис. 1: А(Х1=0, Х2=0, Х3=0); В(8, 0, 2); С(4, 6, 10); А (0, 0, 12); В'(8, 0, 15); С'(4, 6, 18). Плотность тела в вершине А равна 2 г/см3, в вершине С' - 3 г/см3. Координаты вершин призмы заданы в км.
Профиль I совпадает сребром АВ. При этом точка 1 расположена на середине ребра; точка 2 расположена в вершине В; точка 3 - на продолжении ребра в 1 м от вершины В; точки 4, 5 -на продолжении ребра, на удалении от тела.
Профиль II проходит через тело. При этом точки 1, 2 расположены вне тела; точка 3 - на грани АСС'А'; точки 4, 5 - внутри тела; точки 2 и 4 - в 1 м от грани АСС'А' вне и внутри тела, соответственно.
Профиль III проходит над телом параллельно оси ОХх (Х3=0, Х2=0,3; первая точка Х=4, последняя -Х=2004; Ь - расстояние до тела, В - диаметр тела).
В табл. 1, 2 приведены результаты расчетов, выполненных по аналитическим формулам (2)-(5) в сопоставлении с аналогичными расчетами, полученными методом численного интегрирования выражения (1) для Ух3, а также расчетами, приведенными в работе [7]. Численное интегрирование осуществлялось методом Ромберга в системе МаШсаё.
Результаты расчетов тестового примера показывает, что алгоритм, реализующий аналитическое решение задачи (1) устойчив во всех характерных областях: на ребре, в вершине, в плоскости грани, на продолжении ребер и граней, вблизи особых точек, внутри тела.
Время счета алгоритмом (2)-(5) одной точки для одного тела составляет величину порядка 8,7-10-6 с. Расчеты методом численного интегрирования в МаШсаё выполняются на 4 порядка медленнее, не дают результатов в особых точках и в точках, близким к особым. Расчеты по аналитическим формулам в работе [7] проводились на ЭВМ ЕС-1022.
Таблица 2. Результаты вычисления вертикальной составляющей гравитационного потенциала в точках профиля III, расположенных на разных расстояниях от тела
Расстояние от тела Расчет по формулам (2)-(5), 10-5м/с2 Расчет методом численного интегрирования (в МайИса^, 10-5м/с2 Расчет из работы [7], 10-5м/с2
І, км 1/0
0 0 5,98105801-101 5,98105801-Ю1 5,98105801-Ю1
4 0,5 4,10351050-Ю1 4,10351050-Ю1 4,10351050-Ю1
10 1,25 1,48639100-Ю1 1,48639100-Ю1 1,48639100-Ю1
20 2,5 3,69388528-100 3,69388528-100 3,69388528-100
50 6,25 3,21000192-10-1 3,21000192-10-1 3,21000192-10-1
100 12,5 4,21741987-10-2 4,21741987-10-2 4,21741987-10-2
200 25 5,33427875-10-3 5,33427875-10-3 5,33427875-10-2
1000 125 4,27944428-10-5 4,27944428-10-5 4,27945062-10-5
2000 250 5,34895851-10-6 5,34895847-10-6 5,34935387-10-6
Время расчета,с 8,2-10-5 8,0-10-1 13,2
Постановка задачи динамики
неоднородной сильно вязкой среды
Образование и развитие геологических структур во многом обязаны гравитационному тектоге-незу, т. е. движениям, обусловленным аномальными плотностями этих структур на фоне механически равновесного распределения плотности. Плот-ностные неоднородности тектоносферы, наряду с другими тектоническими факторами, вносят свой вклад в напряженно-деформированное состояние среды.
Е.В. Артюшковым [12] на основе анализа послеледниковых изостатических поднятий в Фенноскан-дии и других геолого-геофизических материалов, было показано, что в условиях чрезвычайно медленно протекающих во времени геологических процессов (М03...107лет) и больших размеров геологических тел (Х~103...106м) тектоносферу Земли можно считать сильно вязкой несжимаемой средой с вязкостью т~1019...1024Па-с. Учитывая это, движения среды в тектоносфере могут быть описаны уравнением Навье-Стокса в приближении малых чисел Рейнольдса и уравнением неразрывности, которые в этих условиях принимают следующий вид:
-Ур(х) + Т2 и (х) + Н(х) = 0, V и(х) = 0. (5)
Здесь и(х) - вектор скорости смещения среды в точке х; р(х) - давление в точке х, И(х)=р(х)-£; g=g•i3 - вектор силы тяжести; g=\g\=9,8 м-с-2; ц -орт оси ОХ3; р(х) и т - соответственно плотность и коэффициент динамической вязкости среды.
Исследования, проведённые В.Ю. Косыгиным
[13], Л.А. Масловым и О.С. Комовой [14], показывают, что при моделировании глубинных геодина-мических процессов необходимо учитывать неоднородность строения Земли как по плотностным, так и по реологическим параметрам.
Для того, чтобы учесть эти факторы, моделируемый объем ДсЕ3 будем полагать состоящим из N
N
подобластей (вертикальных призм) Дп: О = и О >
П=1
имеющих общие поверхности контакта йФк.
Пусть т , - значения вязкости среды в каждой из подобластей Д„. Тогда, систему уравнений движения среды можно записать в виде:
^р(х) + тУ2и (х) + Н(х) = 0, V и(х) = 0, хе О , (6)
где Н(х)=рп^, хеДп, и(х) - значение вектора скорости среды; р(х) и рп - соответственно давление и плотность среды, взятые относительно гидростатических значений.
Предполагается, что поверхность дБ моделируемого объема Д свободна от нагрузки, т. е.:
Тх(т,и(х),р(х),п(х)) = 0, х едО. (7)
Здесь символом Тх(т,и(х),р(х),п(х)) обозначен вектор поверхностных сил (напряжений), действующий на бесконечно малый элемент поверхности дБ в точке поверхности х с нормалью п(х):
TЛт, и(хX Р(х), п(х)) = \\*л (т и(х), Р(х))||3х3 •п(х),
||т;1(т,и(х),р(х))||3х3 - матрица, элементами которой являются компоненты напряжений, определяемые соотношениями:
(т и( х) р(х)) =
= т (дмг (х) / дх] + ды} (х) / дх) - 8у- р( х),
где 5ц - символ Кронекера.
На контактах смежных областей Дп и Дк заданы условия непрерывности значений компонент векторов скорости и напряжений:
и(х-0) = и(х + 0), х е 5^, 5^ = дО>п пдД, (8)
Тх (тк, и (х - 0), р( х - 0), п( х)) =
= Тх (лк, и( х + 0), р( х + 0), п( х)). (9)
Необходимо найти распределение скорости смещения среды и(х) и давления р(х), удовлетворяющие системе уравнений (6) при заданных граничных (7) и контактных (8), (9) условиях.
Решение задачи геодинамики
Используя методы теории потенциала [15, 16], решение системы (6)-(9) может быть записано в виде:
и( х) = и1 (х) + и2 (х) + из( х), р( х) = Р( х) + -р (х) + Рз( х), (10)
где щх) и Рг(х) - объёмные потенциалы вида:
и,(х) =т_1(х)£ /О(|,х)• Н(4ЩУ,
п=1 в
Р(х) = £ /Р(4, X) • над;V.
(11)
п=1 В
и2(х) и Р2(х) - потенциалы простого слоя вида: и2(х) — 2г,-\X) ]■ О(х,4) •ъвЩБ,
дВ
Р2(х) = 2 |Р(х,4)<р0(4)(12)
дВ
из(х) и Р3(х) - потенциалы двойного слоя вида:
N
и3(х) = - X (Лп-Лк)Л~1(х)х
п=1,к=п+1
Те(О(4,х),Р(4,х),пп(4))^(4)4 б
(13)
Р3(х) = - X 2(Лп-Лк) х
п=1, к=п+1
/ Егаё х(р(4, х) •п(4)) • ф(4) ^
(14)
л(х)=л„, хєБ„, 0(4,х) - матрица, элементами которой являются величины
^ (4, х) =
= -(8ж)-1 д2В(4,х)/д4 /д^ + д,(4ж)-1 В\4, х), (15)
Б
Б
Р(4,х)=(р'(4,х), р2(£х), р3(4,х)) - вектор, компоненты которого имеют вид:
р (4, х) = -(4я)-1дЯ-1(4, х)/д4,. (16)
Значения функций р(4), 4едД и р(4), 4е8тш определяются из решения системы интегральных уравнений Рисса-Шаудера:
Р0(х) + 2 |Тх(О(х 4), Р(х,4Х п(4)) • %(4)^ =
дО
= /,(х), х е дУ,
N
р (х) + 2 (Тп -т)(Т, +Тр)-1 х
п=1,к=п+1
х| ^ О(4, х), Р(4, х), пп(4))р(4)^ 5 = /(х), (17)
где
/о (х) = “Гх (^1 (x) P1 (х)п(х)) --Гх (^(х), P2(х), И(х)),
Х е Smp — Sконт ;
/1(х) = 2 JG(х,4)(4)^S-
5D
-2(Лт +Лр)£ JG(4,х)• H(4)^4V,
п=1 о
N
5 = и дО \ дО.
конт п
п=1
Решение системы уравнений (6)-(9) в виде функционально-аналитических соотношений (10)-(14) получено Ю.В. Пятаковым [17].
Реализация решения связана с необходимостью вычисления интегралов в правых частях (11)-(14). Структура этих интегральных операторов во многом аналогична структуре интегральных операторов, рассматриваемых в задачах гравитационного потенциала. Поэтому для вычисления правых частей в соотношениях (11)-(14) целесообразно использовать технику решения прямых задач гравиметрии, рассмотренную в первой части статьи. Покажем это на примере вычисления составляющих решения и1(х), Р1(х), определенных соотношениями (11).
Пусть Дп - вертикальная треугольная призма (рис. 2), имеющая избыточную плотность р=рп. Тогда вектор скорости
и( х) = т 1 / О (4, х) • Н (4)^У (18)
Оп
и функция давления
р (х) = / Р (4, х) •Н (4) ¿У (19)
Оп
будут представлять собой решение системы уравнений
-Ур( х) + Т72 и (х) + Н (х) = 0, (20)
V и (х) = 0, (21)
описывающей движение однородной вязкой среды Т=еош1;, обусловленное влиянием тела Дп с избыточной плотностью рп. В соотношениях (18)—(20)
H (х) =
Рп • g, х е Dn
0,
(22)
Подставляя в соотношения (18,19) выражения
(15), (16), (22) для элементов G(4,x), P(4,x) и H(x), получим аналитические выражения для компонент вектора скорости смещения среды u(x) и функции давления p(x):
um (х) = ~Рпё(8^)Z[d^3(m Vi (х) +ei (Х)]7з
'3%
т = 1,2;
и3 (х) =
= -Png (8^Л)-1
Y^vi(х) + ri(х)]^
,=1
+2V (х)
(')
р( х)=-Png (4^)-1 Z Y3(;3v(х).
(23)
(24)
(25)
где
е(х)=Z
а
j=1
Y Y_) / 2
/ 2,3 ' ■‘■J
Г(х), у;(х), г;(х), ау, 4, Р 2,2 - те же, что и в вы-
ражениях (2)-(5) для расчета составляющей Уц(х) гравитационного потенциала.
Из соотношений (23)—(25) видно, что все составляющие аналитических выражений для компонент геодинамического поля и(х) и р(х) полностью определены соответствующими аналитическими составляющими решения прямой задачи гравиметрии.
Тестирование алгоритма решения
прямой задачи геодинамики
В качестве тестового примера рассматривается аномальное тело (рис. 2, а) - горизонтальная призма, моделирующая линзу разуплотненной астеносферы в акватории Охотской глубоководной котловины [18]. Вязкость т вмещающей среды (астеносферы) и самого тела постоянна и равна 1019Па-с, аномальная плотность тела р=-0,02 г/см3. Горизонтальное (верхнее) основание призмы расположено на глубине 100 км перпендикулярно направлению вектора ускорения силы тяжести g. Размеры призмы по оси ОХ1 составляют 500 км, по оси ОХ3 - 250 км, а по оси ОХ2 - 1000 км. Исходное тело разбивается системой типовых аппроксимирующих тел (рис. 1) - вертикальных призм так, как это показано на рис. 2, б.
Проверим, удовлетворяет ли алгоритм, реализующий решение (23)—(25), системе уравнений
(20)-(22)?
S
nk
i =1
Рис. 2. К решению прямой задачи геодинамики: а) расположение тела аномальной плотности; б) аппроксимация тела системой вертикальных призм
Для этого рассмотрим дифференциальный оператор А(т,и(х),р(х)):
А(т, и (х), р(х)) =-р(х) + г(Ч1и (х).
Построим аппроксимацию А(т,и(х),р(х)) разностным оператором А (т,и(х),р(х),А), компоненты которого имеют вид:
Ак = -Л-1[р(х + Л /2- ¡к)- р(х-Л /2- к )]+
3
+ТЛ-2[Еик (х + Л/2-гу.)+ ик (х-Л/2- )-6ц (х)],
1=1
где ц - к-й орт системы координат (к=1,2,3); ик и р - компоненты вектора скорости и давления, вычисленные по формулам (23)-(25); Л - приращение аргумента.
Тогда, вследствие непрерывности и(х), р(х), представленных соотношениями (23)-(25), при Л—0 должно выполняться предельное соотношение вида: А (т,и(х),р(х),Л)—А(т,и(х),р(х)). В этом случае из уравнения (20) следует, что при Л—0 должно выполняться предельное соотношение вида:
А(т, и(х), р(х), Л) + Н(х) — 0. (26)
Проверку соотношения (26) выполним расчетом в точках, расположенных вне и внутри тела (рис. 2, а). Результаты представлены в табл. 3.
Аналогично, заменим оператор Vu(x) в (21) его конечно-разностной аппроксимацией:
ё (и (х), Л) =
=д-1 е (цк(х+Л/2-к )-ик(х-Л/2- ¡к х>.
к=1
Тогда при Л——0 в соответствии с уравнением
(21) должно выполняться предельное соотношение
Таблица 3. Результаты расчета Л(т,и(х),р(х),Л) +Н(х) в точках, расположенных вне и внутри тела
Рас- чет- ные точки Координаты расчетной точки, км Компо- ненты Д(т,и(х), р(х),Л) + Н(х) Численные значения компонент Д(т,и(х),р(х),Л) +Н(х), при разных величинах приращения Л
*1 *2 *3 Л=100 км Л=1 км Л=0,01 км
Вне тела 0 0 100 А -1,0-100 1,9-10-7 -1,0-10-10
А3 1,0-100 1,9-10-7 -3,2-10-11
100 0 500 А 4,6-10-1 4,4-10-8 4,7-10-11
А3 -6,4-10-4 -5,7-10-8 -5,8-10-11
450 0 50 А 0,0 0,0 0,0
А3 -8,5-10-1 2,8-10-7 1,8-10-10
Вну- три тела 300 0 250 А 9,9-10-2 9,3-10-7 1,8-10-10
А3 6,8-10-2 8,1-10-7 0,0
450 0 325 А 0,0 0,0 0,0
А3 2,1-10-2 2,4-10-6 0,0
450 0 350 А 0,0 0,0 0,0
А3 2,0-10-1 1,3-10-5 0,0
ё (и (х), Л) — 0.
(27)
Как следует из результатов, приведенных в табл. 3 и 4, решение, представленное соотношениями (23)-(25), с высокой точностью удовлетворяют исходной системе уравнений (20)-(22).
Тестирование решения, результаты которого приведены в табл. 3 и 4, позволяет выявить возможные ошибки (неточности) в определении по-динтегральных функций, включающих разностные соотношения, в которых фигурируют операции вычисления близких по значению величин.
В качестве дополнительного тестового примера приводим результаты расчетов (табл. 5) компонент вектора и(х) и функциир(х) по формулам (23)—(25) и методом численного интегрирования выражений (18,19), выполненных в точке, внешней по отношению к телу (рис. 2, а). Координаты расчетной точки (в км): Х;=50, х^50, х'3=50.
Рис. 3. К расчету геодинамического поля тектоносферы Охотоморской глубоководной котловины: а) распределение векторов мгновенной скорости смещения u(x) среды, (1ита*1=8,7м/год); б) распределение избыточного давления р(х) в МПа. Пояснения в тексте
Таблица 4. Результаты расчета 6(и(х),Л)хц
Рас- четные точки Координаты расчетной точки,км Численные значения, d(u(x),Д)хr|, при разных значениях приращения Л
*1 *2 *3 Л=100 км Л=1 км Л=0,01 км Л=0,0001 км
Вне тела 0 0 100 -5,0-ТО1 2,2-10-5 2,2-10-9 4,9-10-11
100 0 500 2,3-Ю1 2,0-10-5 2,0-109 3,0-10-11
450 0 50 3,5-Ю1 -1,9-10-5 19-10-9 -11-10-12
Внутри тела 300 0 250 -2,3-100 -11-10-4 -1,1-10-8 -2,8-10-11
450 0 325 1,2-100 1,4-10-4 1,4-10-8 1,3-10-12
450 0 350 1,8-100 5,1-10-4 5,1-10-8 2,2-10-11
Таблица 5. Результаты тестового расчета геодинамических параметров тектоносферы - компонент вектора скорости смещения и(*) и давления р(*). Точка расчета внешняя и удаленная по отношению к аномальному телу
Компоненты и(х), р(х) Расчет по аналитическим формулам (23)-(25) Расчет методом численного интегрирования (в Майлса^
и, м/год -9,60420146-10-1 -9,60420146-10-1
и2, м/год 5,64268622-10-2 5,64268622-10-2
и3, м/год -3,56608779-100 -3,56608779-10°
р, МПа 1,73527102-100 1,73527102-10°
Время расчета, с 9,1-10-5 1,1-10-2
На рис. 3 представлены результаты расчетов поля мгновенных скоростей смещения среды и(х) и давления р(х), обусловленные влиянием аномального тела, взятого в качестве тестового примера (рис. 2, а). Расчеты проводились в плоскости симметрии аномального тела (х2=0). Количество расчетных точек 176, время расчета 1,6-10-2 с.
Полученного распределения скорости смещения среды не противоречит следствиям геодина-мической концепции [12] в зоне субдукции. Ко-
нечно, полученные результаты не вполне корректны, т. к. для примера рассчитаны составляющие только объемных потенциалов (11). Кроме того, нужно учитывать, что над астеносферой расположен 100 км слой литосферы, вязкость которого «1021Па-с. Этой слой естественно и по определению (18) на два порядка «гасит» значения скорости.
Используя различные композиции из типовых аппроксимирующих тел можно со сколь угодной точностью строить модели различных участков тектоносферы Земли и рассчитывать их геодина-мическое состояние.
Выводы
1. На основе типовой методики получено новое аналитическое решение прямой задачи гравиметрии для типового аппроксимирующего элемента - вертикальной треугольной призмы с линейным законом изменения плотности от глубины. На тестовых примерах показаны устойчивость, точность и быстродействие решения в точках, внешних и внутренних по отношению к среде.
2. Решение прямой задачи геодинамики представлено совокупностью потенциалов объемных, простого и двойного слоя. Численная реализация задачи выполнена с использованием типового решения прямых задач гравиметрии. На тестовых примерах продемонстрированы устойчивость, точность и быстродействие алгоритма решения прямой задачи геодинамики.
3. Полученные алгоритмы решения прямых задач могут быть использованы как универсальные для обработки данных наземной, шахтной и скважинной гравиразведки, а также для построения геодинамических моделей тектонос-феры.
СПИСОК ЛИТЕРАТУРЫ
1. Старостенко В.И. Гравитационное поле однородных п-уголь-ных пластин и порождаемых ими призм: обзор // Физика Земли. - 1998. - № 3. - С. 37-53.
2. Страхов В.Н., Лапина М.И. Прямая и обратная задачи гравиметрии и магнитометрии для произвольных однородных многогранников // Теория и практика интерпретации гравитационных и магнитных полей в СССР. - Киев: Наукова думка, 1983. - С. 3-87.
3. Страхов В.Н., Лапина М.Н. Прямые задачи гравиметрии и магнитометрии для однородных многогранников // Геофизический журнал. - 1986. - Т. 8. - № 6. - С. 20-31.
4. Старостенко В.И., Легостаева О.В. Прямая задача гравиметрии для неоднородной произвольно усеченной вертикальной прямоугольной призмы // Физика Земли. - 1998. - № 12. -С. 31-44.
5. Исаев В.И., Гуленок Р.Ю., Лобова ГА. Интерпретация данных высокоточной гравиразведки. Трехмерность объектов // Известия Томского политехнического университета. - 2012. -Т. 320. - №1. - С. 98-104.
6. Пятаков Ю.В., Исаев В.И. Методы решения прямых задач гравиметрии // Известия Томского политехнического университета. - 2012. - Т. 320. - № 1. - С. 105-110.
7. Исаев В.И., Пятаков Ю.В. Решение прямой задачи гравиметрии для трёхмерных блоково-градиентно-слоистых сред // Геофизический журнал. - 1990. - Т. 12. - № 3. - С. 72-79.
8. Голиздра ГЯ. Основные методы решения прямой задачи гравиразведки на ЭВМ // Региональная, разведочная и промысловая геофизика. - М.: ВИЭМС, 1977. - 98 с.
9. Кравцов Г.Г. Поле притяжения многогранников переменной плотности // Записки Ленинградского горного института. -1978. - Вып. 66. - С. 8-17.
10. Балк П.И., Балк ТВ., Носырев В.И. Об аналитическом решении трехмерной прямой задачи гравиразведки в случае переменной плотности возмущающих масс // Известия вузов. Сер. Геология и разведка. - 1976. - № 4. - С. 121-129.
11. Программы по математическому обеспечению обработки и интерпретации геолого-геофизических материалов на ЭВМ / под ред. В.Н. Яковлева. - Л.: Изд-во НПО «Рудгеофизика», 1982. - 354 с.
12. Артюшков Е.В. Геодинамика. - М.: Наука, 1979. - 328 с.
13. Косыгин В.Ю. Гравитационное поле и плотностные модели тектоносферы Северо-Запада Тихого океана. - Владивосток: ДВО АН СССР, 1991. - 201 с.
14. Маслов Л.А., Комова О.С. Численное моделирование глубинных геодинамических процессов в активных окраинах // Физика Земли. - 1990. - № 3. - С. 53-60.
15. Михлин С.Г. Курс математической физики. - СПб.: Изд-во «Лань», 2002. - 576 с.
16. Трёхмерные задачи математической теории упругости и термоупругости / под ред. В.Д. Купрадзе. - М.: Наука, 1976. - 664 с.
17. Пятаков Ю.В. Трехмерная задача динамики сильно вязкой несжимаемой многокомпонентной гравитирующей среды. I. Постановка и алгоритм решения задачи // Геофизический журнал. - 2005. - Т. 27. - № 3. - С. 387-392.
18. Красный М.Л., Косыгин В.Ю., Исаев В.И. Оптимальная плот-ностная модель тектоносферы вдоль геотраверса о. Сахалин -о. Итуруп - Тихий океан // Тихоокеанская геология. - 1985. -Т. 4. - № 6. - С. 36-48.
Поступила 27.04.2012 г.
