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

Используя эту подстановку, получим:u = w − u + q , u = w − 2w u + q + u + q ,2u=w−2222222Эйлера:w2 + q2,u +q =2w22w2 + q22w ⋅ 2w − 2( w 2 − q 2 )w2 + q2 w2 − q2=dw, du =dw .=2w2w4w 22w 261Поскольку теперь интегрирование происходит по переменной w, тоизменятся и пределы интегрирования: нижний предел – (η1ν − y ν ) + R1 ,верхний – (η ν2 − yν ) + R2 .Дальнейшие преобразования будут следующими:ην2 − yν∫η1ν − yνdup + u2 + q 2=(ην − yν )+ R22∫(ην − yν )+ R11w2 + q2(ην2 − yν )+ R22w( w 2 + q 2 )2w 2dw =∫ 2w 2 ( w 2 + 2 pw + q 2 ) dw =w2 + q2(η1ν − yν )+ R1p+2w(ην − yν )+ R(ην − yν )+ R22⎛1⎞2pw 2 + 2 pw + q 2 − 2 pw⎜⎜ − 2⎟dw =dw ==222 ⎟∫∫(2)2ww+pw+qww+pw+qνννν⎝⎠(η1 − y )+ R1(η1 − y )+ R122(ην − yν )+ R(ην − yν )+ R2211=dw−2pdw = I 21 − 2 pI 22 .∫∫22(η1ν − yν )+ R1 w(η1ν − yν )+ R1 w + 2 pw + q22Первый интеграл легко беретсяI 21 =(ην − yν )+ R(η ν2 − yν ) + R2 .1dw=ln∫(η1ν − yν ) + R1(η1ν − yν )+ R1 w22Рассмотрим второй интеграл:I 22 =(ην − yν )+ R1dw .2wpwq+2+νν(η1 − y )+ R12∫22Для того, чтобы его взять введем следующую подстановку z = w + p .Тогдаw 2 + 2 pw + q 2 = ( z − p ) 2 + 2 p( z − p ) + q 2 = z 2 − p 2 + q 2 .С учетом того, что ζ − z = p , (ξ ν − xν ) 2 + (ζ − z ) 2 = q 2 , последнеевыражение приобретет вид:w 2 + 2 pw + q 2 = z 2 − p 2 + q 2 = z 2 + (ξ ν − xν ) 2 ,и, соответственно,62(ην − yν )+ R(ην − yν )+ R + p22dz1dw=∫ν w 2 + 2 pw + q 2∫ z 2 + (ξ ν − xν ) 2 .ννν(η1 − y )+ R1(η1 − y )+ R1 + p22Сделаем еще одну подстановку: v =(ην2ν(ην − yν )+ R + p2)− y + R2 + p∫dz(ην − yν )+ R + p z + (ξ − x )1z, тогда(ξ − zν )νν2ν2=1ξ ν − xν21dv.ν∫ξ − x (η1ν − yν )+ R1 + p 1 + v 2νξ ν − xνПолученный интеграл представляет собой разность арктангенсов, и, впринципе, на этом можно было бы остановиться.

Однако мы продолжимпреобразования. Это связано с тем, значения арктангенсов зависят отрасстояний R1 и R2. При их больших величинах значения арктангенсовбудут близки к (π/2), в то же время значения самого интеграла сувеличением R1 и R2 должно убывать. Это означает, что мы можемпотерять точность при вычислении значений поля по этой формуле.Поэтому преобразуем полученный интеграл, введя очереднуюподстановку:1v= ,v′1dv = −dv ′ ,( v ′) 2−1dv ′(v ′) 2⎛1⎞1+ ⎜ ⎟⎝ v′ ⎠(ην − yν )+ R + p2νν2ξ −x1dv1= ννν∫2ξ − x (η1ν − yν )+ R1 + p 1 + v ξ − xνννξ −x1= νξ − xν(ηξ ν − xνν1(η2⎛1 ⎞ ( v ′) 2dv ′′⎟=−;dv= ⎜⎜ −2 ⎟21 + ( v ′) 2⎝ ( v ′) ⎠ 1 + (v ′))− yν + R1 + p∫ξ ν − xνν2)dv ′=1 + ( v ′) 2− yν + R2 + p⎛⎞ξ ν − xνξ ν − xν⎜⎜ arctg ν⎟−arctg(η1 − yν ) + R1 + p(η ν2 − yν ) + R2 + p ⎟⎠ .⎝Введя обозначения:w1ν = (η1ν − yν ) + R1 + p ,wν2 = (η ν2 − yν ) + R2 + p , где p = ζ − z ,63и записав разность арктангенсов в виде одного арктангенса, последнеевыражение можно представить в виде:⎞ξ ν − xνξ ν − xν1 ⎛⎜ arctg ν⎟− arctg ννν ⎜νν(η1 − y ) + R1 + p(η 2 − y ) + R2 + p ⎟⎠ =ξ −x ⎝(ξ ν − xν )( wν2 − w1ν )1= νarctg ν ν.ξ − xνw 2 w1 + (ξ ν − xν ) 2При этом мы воспользовались соотношениемarctgX − arctgY = arctgЭтосправедливоξ − xν(η1ν − yν ) + R1 + p именьше 1.νX −Y, при условии, что XY > −1.1 + XYдлянашего случая, поскольку выраженияξ − xν(ην2 − yν ) + R2 + p всегда по абсолютной величинеνНапомним, с чего начинались наши преобразования.

Мыпредставили вертикальную составляющую поля силы притяжениягоризонтальной пластины в следующем виде:1 N ν(g z ( M 0 ) = Gδ пξ − xν )(I1 − I 2 ) .∑ζ − z ν =1Первый интеграл, как мы уже отмечали, равенην2∫rη1ν1dηMM 0ν(η= ln(ην2ν1− y ν ) + R2.− yν ) + R1Второй интеграл после преобразований оказался равен разности двухинтегралов:I 2 = I 21 − 2 pI 22 ,первый из которых равенI 21(η= ln(ην2ν1− yν ) + R2,− yν ) + R164а второй –(ξ ν − xν )( wν2 − w1ν )1I 22 = νarctg ν ν.ξ − xνw 2 w1 + (ξ ν − xν ) 2Поскольку интегралы I1 и I21 равны (I1 = I21), то выражение длявертикальной составляющей приобретает вид:1 N ν(g z ( M 0 ) = Gδ пξ − xν )2 pI 22 =∑ζ − z ν =1(ξ ν − xν )( wν2 − w1ν ) ⎞1 N ν1ν ⎛= Gδ п∑ (ξ − x )⎜⎜ 2 ζ − z ξ ν − xν arctg wν wν + (ξ ν − xν ) 2 ⎟⎟ =ζ − z ν =1⎝⎠2 1(ξ ν − xν )( wν2 − w1ν )= 2Gδ п sign(ζ − z )∑ arctg ν ν.w 2 w1 + (ξ ν − xν ) 2ν =1N13. Теперь можно снова вернуться к выражению для потенциалапластины.

Как мы уже отмечали, это выражение может быть представленов следующем виде:NVq ( M 0 ) = δ п ∑ (ξ ν − x ν )V лν ( M 0 ) − δ п (ζ − z )ν =1∂V q ( M 0 )∂z,гдеR1ν + Rν2 + LνV л ( M 0 ) = G ln νR1 + Rν2 − Lνν– потенциал отрезка, совпадающего с ν–ой стороной пластины, сединичной линейной плотностью, R1ν и Rν2 – расстояния от концов отрезкадо точки наблюдения, Lν – его длина;∂ Vq ( M 0 )∂z(ξ ν − xν )( wν2 − w1ν )= 2G sign(ζ − z )∑ arctg ν νw 2 w1 + (ξ ν − xν ) 2ν =1N– вертикальная составляющая притяжения пластины с единичнойповерхностной плотностью;w1ν = (η1ν − yν ) + R1ν + ζ − z , wν2 = (η ν2 − yν ) + Rν2 + ζ − z .65Значения (xν,yν,z) соответствуют координатам расчетной точки M0, азначения (ξ ν ,η1ν ,ζ ) и (ξ ν ,η ν2 ,ζ ) – координатам ν–ой сторонымногоугольника в системе координат, связанной с этой стороной, при этомось oXν направлена ортогонально ν–ой стороне и лежит в плоскостипластины, а ось oYν параллельна этой стороне.

Система координат oXνYνZобразует правую тройку.14. Подведем некоторый итог. Нами были получены выражения дляпотенциала силы притяжения и его производных для однородногомногогранника:V (M0 ) =δQ∑ (ζ2 1qq=− z q )Vq ( M 0 ) ;Qrrg ( M 0 ) = −δ ∑1nq Vq ( M 0 ) .q =1Как видно из этих выражений, поле силы притяжения такогомногогранника определяется через потенциалы его граней, впредположении, что эти грани представляют собой многоугольныепластины с постоянной единичной плотностью.В свою очередь потенциал пластины определяется через потенциалыматериальных отрезков совпадающих с ребрами пластины и компонентуполя силы тяжести ортогональную к плоскости пластины:NVq ( M 0 ) = δ п ∑ (ξ qν − x qν )V лν ( M 0 ) − δ п (ζ q − z q )ν =1∂ Vq ( M 0 )∂z q.Таким образом, алгоритм расчета поля от многогранника будетследующим.

Для каждой q–ой грани многоугольника вводятся новыесистемы координат так, что ось oZq совпадает с внешней нормалью к этойстороне, а оси oXqν и oYqν определяются сторонами грани. Вопрос о том,как вводится такая система координат, мы рассматривали в предыдущейлекции. В этой новой системе координат рассчитываются координатырасчетной точки M0(xqν,yqν,zq) и координаты начала (ξ qν ,η1qν ,ζ q ) и конца(ξ qν ,η 2qν ,ζ q ) ν–ой стороны этой многоугольной пластины. Далеевычисляются значения потенциалов отрезков, совпадающих со сторонамипластины, и вертикальная составляющая притяжения самой пластины вэтой системе координат, после чего определяется потенциал самойrпластины в точке M0.

Зная направление нормали nq к q–ой грани,определяются компоненты поля силы притяжения, связанные с q–ойгранью многогранника, уже в исходной, основной, системе координат66oXYZ. Для того чтобы получить эффект от всего многогранниканеобходимо будет рассчитать эффекты от всех граней и ихпросуммировать.15. Перейдем к вопросу вычисления аномальногомагнитного поляrмногогранника с постоянной намагниченностью I . Как нами ранее былопоказано, решение этой задачи может быть основано на двух подходах. Воснове одного из них лежит соотношение Пуассона о связи потенциаловгравитационного и магнитного полей, второй – на представлении поля ввиде эффектов поверхностных фиктивных “магнитных зарядов”.Так, соотношение Пуассона для компонент аномального магнитногополя имеет вид:⎡V xx V xy V xz ⎤ ⎡ I x ⎤⎡X ⎤µ⎥⎢ ⎥⎢ Y ⎥ = 0 ⎢VVVyxyyyz⎥ ⎢I y ⎥ ,⎢ ⎥ 4π ⎢⎢Vzx Vzy Vzz ⎥ ⎢⎣ I z ⎥⎦⎢⎣ Z ⎥⎦⎦⎣где Ix, Iy, Iz – компоненты вектора намагниченности, Vij – вторыепроизводные потенциала V, определяемого следующим образом:V (M0 ) = ∫D1rMM 0dv .С учетом того, что компоненты аномального гравитационного поля,создаваемого многогранником, выражаются через потенциалы его гранейVq:Q∂V ( M 0 )= − ∑ cos( z q , ri )Vq ( M 0 ) ,∂riq =1то вторые производные будут представляться через компоненты поля,создаваемого пластинами, совпадающими с гранями многогранника иимеющими единичную поверхностную плотность:Q∂V q ( M 0 )∂ 2V ( M 0 ).= − ∑ cos( z q , ri )∂r j∂ri ∂r jq =1В этих формулах ri и rj могут принимать значения x, y, z; Q – число гранейв многограннике, направление оси zq совпадает с направлением внешнейнормали к q–ой грани многогранника.6716.

Второй подход основан на том, что магнитное поле однороднонамагниченного объема можно представить как поле некоторыхфиктивных “магнитных зарядов”, сосредоточенных на поверхности этогообъема. В частности для скалярного магнитного потенциала нами былоранее получено следующее соотношение:r rI ⋅ 1n∫S rMM dS .01 δ s (M )1dS =U(M0 ) =∫4π S rMM 04πПрименительно к многограннику это соотношение приобретет вид:1U(M0 ) =4πr rI ⋅ 1n1dS=∫S rMM4π0Q∑1 ∫q = Sqr rI ⋅ 1nqrMM 01dS =4πδqQ∑1 ∫ rq = Sq MM 01dS =4πQ∑1 δ V ( M 0 )q=qq,r rгде δ q = I ⋅ 1nq – поверхностная плотность “магнитных зарядов” q–ойграни, Vq – потенциал притяжения q–ой грани.Напряженность аномального магнитного поля определяется черезградиент скалярного потенциала:rH ( M 0 ) = − gradU ( M 0 ) .Для компонент поля создаваемого многогранником, и с учетом того, чтонас интересует не значения напряженности, а значения индукции,получим:µ∂U ( M 0 )=− 0X ( M 0 ) = −µ 0∂x4πY ( M 0 ) = −µ 0µ∂U ( M 0 )=− 0∂y4πµ∂U ( M 0 )=− 0Z ( M 0 ) = −µ 0∂z4πQ∑δ qq =1q=qq=∂y∂ Vq ( M 0 )Q∑1 δ∂x∂ Vq ( M 0 )Q∑1 δ∂ Vq ( M 0 )q∂z,,.17.

Таким образом, как при вычислении компонент аномальногомагнитного поля на основе соотношения Пуассона, так и при вычисленииполя через фиктивные “магнитные заряды”, алгоритм будет один и тот же.Вводится новая система координат, связанная с q–ой граньюмногогранника. Ось oZq этой системы будет совпадать по направлению свнешней нормалью к этой грани. Далее, на основании полученных ранее68соотношений для элементов притяжения горизонтальной многоугольнойпластины, рассчитываются компоненты ее силы притяжения. Если расчетведется через фиктивные “магнитные заряды”, то вычисляется и ихповерхностная плотность на каждой грани.

После чего вычисляетсяэффект этой грани в исходной системе координат. Общий эффект будетопределяться как сумма эффектов от каждой грани в отдельности.18. Сделаем еще одно замечание. Казалось бы, что расчет поля черезфиктивные “магнитные заряды” представляется менее трудоемкойоперацией по сравнению с расчетом поля на основе соотношенияПуассона. Тем не менее, по трудоемкости эти вычисления совпадают. Длятого чтобы это показать представим поверхностную плотность зарядов q–ой грани следующим образом:rrrrrδ q = I ⋅ 1nq = I x cos( x , nq ) + I y cos( y , nq ) + I z cos( z , nq ) .Вторые производные потенциала притяжения многогранника выражаютсяследующим образом:Q∂V q ( M 0 )∂ 2V ( M 0 ).= − ∑ cos( z q , ri )∂r j∂ri ∂r jq =1∂ 2V ( M 0 ) ∂ 2V ( M 0 )∂ 2V ( M 0 )С учетом того, что, выражение дляможно=∂ri ∂r j∂r j ∂ri∂ri ∂r jпредставить в виде:Q∂ Vq ( M 0 )∂ 2V ( M 0 ).= − ∑ cos( z q , r j )∂ri ∂r j∂rq =1iВыпишем для примера X компоненту аномального магнитного поля,получаемую из соотношения Пуассона:µ 0 ⎛ ∂ 2V∂ 2V∂ 2V ⎞⎟=⎜+ Iy+ IzX=Ix4π ⎜⎝ ∂x 2∂x∂y∂x∂z ⎟⎠=QQQ∂V∂V∂V ⎞µ0 ⎛⎜⎜ − I x ∑ cos( z q , x ) q − I y ∑ cos( z q , y ) q − I z ∑ cos( z q , z ) q ⎟⎟ =∂x∂x∂x ⎠4π ⎝q =1q =1q =1QQ∂Vq∂ Vq∂Vq ⎞µ0 ⎛ Qqq⎜⎜ ∑ I x cos( z , x )⎟==−+ ∑ I y cos( z , y )+ ∑ I z cos( z q , z )∂x q =1∂x q =1∂x ⎟⎠4π ⎝ q =169∂V ⎞µ0 ⎛ Q⎜⎜ ∑ (I x cos( z q , x ) + I y cos( z q , y ) + I z cos( z q , z )) q ⎟⎟ .=−∂x ⎠4π ⎝ q =1rПоскольку ось oZq совпадает по направлению с внешней нормалью nq , товыражение, стоящее под знаком суммы в круглых скобках будет равноповерхностной плотности зарядов q–ой грани:I x cos( z q , x ) + I y cos( z q , y ) + I z cos( z q , z ) = δ q .Таким образом, для компоненты X получаем:µ 0 ⎛ Q ∂ Vq ⎞⎜ ∑δ q⎟,X =−∂x ⎟⎠4π ⎜⎝ q =1что совпадает с полученным ранее выражением при вычислении полячерез поверхностные заряды.Литература.1.