А.А. Булычев, И.В. Лыгин, В.Р. Мелихов - Численные методы решения прямых задач грави- и магниторазведки, страница 9

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

Далее, по полученной им формуле былирассчитаны их эффекты. В дальнейшем формулы силы тяжести отгоризонтальной многоугольной пластины были получены разнымиисследователями, и с ними можно познакомиться в Справчнике геофизика,том Гравиразведка и другой литературе. Однако выражение длягравитационного потенциала пластины видимо впервые было получено53В.Н. Страховым. Им же были получены новые выражения для поля силытяжести и его производных создаваемых такой моделью.6.Рассмотрим модель горизонтальной пластины SN c N сторонами.

Дляэтого введем декартову систему координат, и будем предполагать, чтопластина параллельна плоскости oXY и лежит на глубине ζ. Поскольку мывновь будем использовать прием, связанный с понижением кратностиинтеграла, то необходимо задать направление обхода пластины. Обходпластины задает нам направление нормали, которое определяется поправилу буравчика. Если направление нормали будет совпадать снаправлением оси oZ, то такой обход будем считать положительным.Потенциалпластиныопределяется следующим образом:Vq ( M 0 ) = Gδ п∫rSN1dS ,MM 0где δп – поверхностная плотностьпластины.

Для нас представляет интерес получить выражения как длясамого потенциала, поскольку через него выражаются элементы полеймногогранника, так и для его производных, посколькуrg ( M 0 ) = gradVq ( M 0 ) .7.Для того, чтобы получить эти выражения воспользуемся подходом,предложенным В.Н. Страховым. Им были сконструированы следующиеконструкции:1rMM 0= ∆ ξ ,η (rMM 0 − ζ − z ln( rMM 0 + ζ − z ),(ζ − z ) 2= ζ − z ∆ ξ ,η ln( rMM 0 + ζ − z ) .3rMM08.Рассмотрим вначале, каким образом можно представить компонентыполя притяжения.

Здесь надо рассмотреть два случая: 1–й – как будутпредставляться горизонтальные компоненты поля в нашей системекоординат, и 2–й – представление вертикальной компоненты поля.Поскольку (ζ–z) для нашей модели величина постоянная, то длягоризонтальных компонент поля притяжения можно записать:54⎛⎞1g x 1x + g y 1 y = grad x , yVq ( M 0 ) = grad x , y ⎜ Gδ п ∫dS ⎟ =⎜⎟S N rMM 0⎝⎠= −Gδ п∫SNgrad ξ ,η1rMM 0dS = −Gδ п1∫∂S NN1 rr1n dl = −Gδ п ∑ ∫1nν dl ,rMM 0ν =1 Lν rMM 0rгде 1nν – орт нормали к ν–ой стороне и лежащей в плоскости пластины.rПоскольку для каждой из сторон он имеет свое направление, то орт 1nνможно записать перед знаком интеграла. Кроме того, выражение1G∫dl представляет собой потенциал материальной линии сrLν MM 0единичной линейной плотностью:1V лν ( M 0 ) = G ∫dl .Lν rMM 0В результате для горизонтальных компонент поля получим:Nrrrg x 1x + g y 1 y = −Gδ п ∑1nνν =1∫rLν1MM 0Nrdl = − δ п ∑1nν V лν ( M 0 ) .ν =1Покомпонентная запись будет иметь вид:gx (M0 ) =gz ( M 0 ) =∂ Vq ( M 0 )∂x∂Vq ( M 0 )∂zNr r= −δ п ∑ cos(1nν ,1x )V лν ( M 0 ) ,ν =1Nr r= −δ п ∑ cos(1nν ,1z )V лν ( M 0 ) .ν =1Стоит обратить внимание на схожесть этих формул с формулой для полясилы тяжести многогранника, выражаемой через потенциалы его граней:QQ1rrrg ( M 0 ) = −Gδ ∑1nq ∫dS = −δ ∑1nq Vq ( M 0 ) .q =1q =1Sq rMM 09.Рассмотрим, каким образом можно представить выражение длявертикальной составляющей поля силы притяжения нашей модели.Выпишем это выражение:55gz ( M 0 ) =∂ Vq ( M 0 )∂z= Gδ п∫ζ −zSN3rMM0dS ,и воспользуемся предложенной В.Н.

Страховым подстановкой(ζ − z )2= ζ − z ∆ ξ ,η ln( rMM0 + ζ − z ) .3rMM0Тогдаg z ( M 0 ) = Gδ п= Gδ пζ −zζ −z11(ζ − z ) 2dS = Gδ пζ − z ∆ ξ ,η ln( rMM 0 + ζ − z )dS =∫3ζ − z S N rMM 0ζ − z S∫N∫ ∆ ξ η ln( r,MM 0+ ζ − z )dS = Gδ пSNζ −zζ −z∫∂ ln( rMM 0 + ζ − z )∂n∂S Ndl ,rгде вектор нормали n , как и в предыдущем случае, лежит в плоскостипластины.10. Теперь вновь вернемся к потенциалу пластины и воспользуемсяподстановкой1= ∆ ξ ,η (rMM 0 − ζ − z ln( rMM 0 + ζ − z ) .rMM 0С учетом этой подстановки потенциал пластины можно представить ввиде:Vq ( M 0 ) = Gδ п∫rSN1MM 0dS = Gδ п ∫ ∆ ξ ,η (rMM 0 − ζ − z ln( rMM 0 + ζ − z )dS =SN= Gδ п∫∂S N∂rMM 0∂ndl − Gδ п ζ − z∫∂ ln( rMM 0 + ζ − z )∂n∂S Ndl .Сравним второе слагаемое полученного выражения с выражением длявертикальной компоненты поля, создаваемой рассматриваемой модельюпластины:g z ( M 0 ) = Gδ пζ −zζ −z∫∂S N56∂ ln( rMM 0 + ζ − z )∂ndl .Видно, что второе слагаемое можно выразить через компоненту gz поляпритяжения пластины с единичной поверхностной плотностью:Gδ п ζ − z∫∂ ln( rMM 0 + ζ − z )∂n∂S Ndl = δ п (ζ − z ) g z ( M 0 ) = δ п (ζ − z )∂V q ( M 0 )∂z.Рассмотрим первое слагаемое в выражении для потенциалапластины, и представим его в виде суммы интегралов по сторонаммногоугольника:Gδ п∫∂S N∂rMM 0∂nNdl = Gδ п ∑ ∫ν =1 Lν∂rMM 0∂nνdl .Для того, чтобы вычислитьпроизводную по нормали к ν–ойстороневоспользуемсяужеопробованным приемом – введемновуюсистемукоординат,связаннуюсосторонамимногоугольника.

В этой системекоординат ось oZ сохранит свое направление, ось oXν будет совпадать с ν–ой нормалью, ось oYν – совпадать с ν–ой стороной. Во введенной такимобразом системе координат расстояние rMM 0 будет определятьсяследующим образом:rMM 0 = (ξ ν − x ν ) 2 + (η ν − y ν ) 2 + (ζ − z ) 2 ,производная∫Lν∂rMM 0∂nνdl =∫Lν∂∂= ν . Тогда∂nν ∂ξ∂rMM 0∂ξ νdl =ην2(ξ ν − xν ) ν∫ν rMM dη =η10ην2= (ξ ν − xν ) ∫1rdη ν = (ξ ν − xν )V лν ( M 0 ) ,η1ν MM 0где V лν ( M 0 ) – потенциал материального отрезка с единичной линейнойплотностью и совпадающего с ν–ой стороной многоугольника.57В результате всех проведенных преобразований для потенциалагоризонтальной многоугольной пластины можем записать:NVq ( M 0 ) = δ п ∑ (ξ ν − x ν )V лν ( M 0 ) − δ п (ζ − z )ν =1∂V q ( M 0 )∂z.Таким образом, потенциал горизонтальной многоугольной пластиныопределяется через потенциалы материальных отрезков, совпадающих ссторонами многоугольника и через вертикальную компоненту силыпритяжения этой пластины.11.

Вновь вернемся к потенциалу притяжения Vл(M0) материальногоотрезка. Ранее уже было получено его выражение для частного случая, аименно, для случая, когда отрезок совпадает по направлению с одной изкоординатных осей. Так для стержня, расположенного вдоль оси oX этовыражение имеет вид:V л ( M 0 ) = Gσ л ln(ξ 2 − x ) + R2,(ξ 1 − x ) + R1где ξ1 и ξ2 – координаты положения концов отрезка по оси oX, R1 и R2 –расстояния от начала и конца отрезка до точки наблюдения M0.Затем нами было показано, как с помощью введения новойкоординатной системы, связанной с этим отрезком, можно получитьвыражения и осуществить вычисления элементов поля, создаваемогоотрезком при его произвольном расположении относительно основнойсистемы координат.Не составляет особого труда ввести систему координат, связанную сэтим отрезком, провести в ней вычисления, и вернуться к исходнойсистеме.

Таким образом, приведенная формула для потенциала вполнеудобна, но тем не менее, она не совсем хороша с вычислительной точкизрения. Для примера рассмотрим случай, когда отрезок лежит на оси oX. Вэтом случае расстояния R1 и R2 будут равны:R1 = (ξ 1 − x ) 2 + y 2 + z 2 ,R2 = (ξ 2 − x ) 2 + y 2 + z 2 .Пусть точка наблюдения M0расположена на оси oX, в точке скоординатами (x,0,0), при этомx>ξ2. Рассмотрим, чему будетравен числитель дроби, стоящийпод знаком логарифма:58(ξ 2 − x ) + R2 = (ξ 2 − x ) + (ξ 2 − x ) 2 = (ξ 2 − x ) + ξ 2 − x = 0 ,поскольку (ξ 2 − x ) < 0 . Соответственно сама дробь оказывается равнойнулю, и логарифм от такой дроби оказывается неопределенным (логарифмнуля – бесконечно большая отрицательная величина).

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

В.Н. Страховым такаяформула была предложена:V л ( M 0 ) = Gδ л lnR1 + R2 + L,R1 + R2 − Lгде R1 и R2 – расстояния от концов отрезка до точки наблюдения M0, L –длина отрезка.Для того чтобы показать справедливость этой формулы повторимпредложенные им выкладки, при этом будем предполагать, что отрезокрасполагается вдоль оси oX. Тогда,ξ2 – ξ1 = L – длина отрезка;(ξ2 – x) – (ξ1 – x) = L;R22 − R12 = ((ξ 2 − x ) 2 + (η − y ) 2 + (ζ − z ) 2 ) − ((ξ 1 − x ) 2 + (η − y ) 2 + (ζ − z ) 2 ) == (ξ 2 − x ) 2 − (ξ1 − x ) 2 = L((ξ 2 − x ) + (ξ1 − x )) ;R22 − R12((ξ 2 − x ) + (ξ 1 − x )) =;L((ξ 2 − x ) − (ξ 1 − x )) = L ;⎞1 ⎛ R22 − R12ξ 2 − x = ⎜⎜+ L ⎟⎟ ,L2⎝⎠⎞1 ⎛ R22 − R12ξ 1 − x = ⎜⎜− L ⎟⎟ ;L2⎝⎠59⎞1 ⎛ R22 − R12⎜⎜+ L ⎟⎟ + R2L2R22 + 2 R2 L + L2 − R12⎠=σln=V л ( M 0 ) = Gσ л ln ⎝ 2Gл2222−+−RL2RLR⎛⎞1 R2 − R1211⎜⎜− L ⎟⎟ + R12⎝L⎠= Gσ л ln( R2 + L − R1 )( R2 + L + R1 )( R2 + L) 2 − R12Gln=σ=л( R2 − R1 + L)( R2 + R1 − L)R22 − ( R1 − L) 2= Gσ л ln( R2 + R1 + L).( R2 + R1 − L)12.

Получим теперь выражение для вертикальной составляющей поляпритяжения горизонтальной многоугольной пластины. Как ужеотмечалось, обход пластины должен быть выбран таким образом, чтонормаль к пластине совпадает по направлению с осью oZ. Для этого случаянами было получено следующее представление для вертикальнойкомпоненты поля gz:g z ( M 0 ) = Gδ пζ −zζ −z∫∂ ln( rMM 0 + ζ − z )∂S N∂ndl .Представим интеграл по замкнутому контуру суммой интегралов поотрезкам, образующим стороны многоугольника, и вновь для каждой ν–ойстороны введем свою систему координат, так, что ось oXν будетнаправлена по нормали к этой стороне, ось oYν будет совпадать с ν–ойстороной, а ось oZ останется без изменений. Тогда, как мы это уже делалиранее, можем осуществить следующие преобразования:g z ( M 0 ) = Gδ пζ −z= Gδ пζ −zζ −zζ −zN∑∫ν=1 LνN∑∫ν∂ ln( rMM 0 + ζ − z )=1 Lν∂nν∂ ln( rMM 0 + ζ − z )∂ξ νdl =ζ −zdl = Gδ пζ −zν(ξνN η2∑∫ (rν=1 ην1MM 0− xν )dη ν =+ ζ − z )rMM 0νη2ζ −z1 N νν(ξ − x ) ∫= Gδ пdη ν =∑( rMM 0 + ζ − z )rMM 0ζ − z ν =1η1ν60η2⎛ η2 1⎞1 N ν1νν ⎜(= Gδ пξ −x ) ∫dη − ∫dη ν ⎟ =∑⎜ ν rMM⎟+ζ −zζ − z ν =1rη1ν MM 00⎝ η1⎠νν1 N ν(= Gδ пξ − xν )(I1 − I 2 ) .∑ζ − z ν =1Первый интеграл этой суммы нами уже рассматривался:I1 =ην2∫r1dηη1ν MM 0ν(η= ln(ην2ν1− yν ) + R2,− yν ) + R1где, как и раньше, R1 и R2 – расстояния от концов отрезка до точкинаблюдения M0.Рассмотрим 2–ой интеграл:I2 =ην2ην∫rη1ν MM 0211dη ν = ∫dη ν ,2νν 2νν 2+ζ −zη1ν ζ − z + (ξ − x ) + (η − y ) + (ζ − z )и введем обозначения:ζ − z = p,(ξ ν − xν ) 2 + (ζ − z ) 2 = q 2 ,(η ν − y ν ) = u .Тогда этот интеграл приобретет вид:ην2 − yνν∫νη1 − ydup + u2 + q2.Такой интеграл можно взять с помощью подстановкиw = u + u 2 + q 2 .