Почему Земля такая жирная?

Ошибка заключается в том, что вы предполагаете, что распределение плотности «почти сферически симметрично». Это достаточно далеко от сферической симметрии, если вы хотите рассчитать субведущие эффекты первого порядка, такие как экваториальная выпуклость. Если ваша цель состоит в том, чтобы вычислить отклонения уровня моря от сферической симметрии (до первого порядка), непоследовательно пренебрегать такими же большими поправками первого порядка к сферической симметрии с другой стороны — источника гравитации. Другими словами, термин$hg$в вашем потенциале неправильно.

Только представьте, что Земля представляет собой эллипсоид с экваториальной выпуклостью, она не вращается и на поверхности нет воды. Каким был бы потенциал на поверхности или потенциал на фиксированном расстоянии от центра эллипсоида? Вы де-факто предполагали, что в этом случае будет$-GM/R+h(\theta)g$где$R$- фиксированный радиус Земли (сферического распределения вещества) и$R+h(\theta)$- фактическое расстояние зонда от начала координат (центра Земли). Однако этим анзацем вы только признали переменное расстояние зонда от сферически симметричного источника гравитации: вы все еще пренебрегли вкладом выпуклости в несферичность гравитационного поля.

Если включить несферически-симметричную поправку к гравитационному полю Земли,$hg$примерно изменится на$hg-hg/2=hg/2$, и, соответственно, искомая выпуклость$\Delta h$придется удвоить, чтобы компенсировать вращательный потенциал. Эвристическое объяснение фактора$1/2$заключается в том, что истинный потенциал над эллипсоидом зависит от «нечто промежуточного» между расстоянием от центра масс и расстоянием от поверхности. Другими словами, «постоянная потенциальная поверхность» вокруг эллипсоидального источника материи находится «точно посередине» между фактической поверхностью эллипсоида и сферической поверхностью.$R={\rm const}$поверхность.

Я постараюсь добавить более точные формулы для гравитационного поля эллипсоида в обновленной версии этого ответа.

Обновление: гравитационное поле эллипсоида

Я численно проверил, что гравитационное поле эллипсоида имеет именно тот эффект деления пополам, который я набросал выше, используя код Монте-Карло Mathematica, чтобы избежать двойных интегралов, которые можно было бы вычислить аналитически, но пока это меня раздражает.

Я взял миллионы случайных точек внутри вытянутого эллипсоида с «радиусами».$(r_x,r_y,r_z)=(0.9,0.9,1.0)$; обратите внимание, что разница между двумя радиусами$0.1$. Среднее значение$1/r$, обратное расстояние между случайной точкой эллипсоида и выбранной точкой над эллипсоидом равно$0.05=0.1/2$меньше, если выбранная точка находится над экватором, чем если она находится над полюсом, при условии, что расстояние от начала координат одинаково для обеих выбранных точек.

Код:

{xt, yt, zt} = {1.1, 0, 0};

runs = 200000;
totalRinverse = 0;
total = 0;

For[i = 1, i < runs, i++,
 x = RandomReal[]*2 - 1;
 y = RandomReal[]*2 - 1;
 z = RandomReal[]*2 - 1;
 inside = x^2/0.81 + y^2/0.81 + z^2 < 1;
 total = If[inside, total + 1, total];
 totalRinverse = 
  totalRinverse + 
   If[inside, 1/Sqrt[(x - xt)^2 + (y - yt)^2 + (z - zt)^2], 0];
]

res1 = N[total/runs / (4 Pi/3/8)]
res2 = N[totalRinverse/runs / (4 Pi/3/8)]
res2/res1

Описание

Используйте приведенный выше код Mathematica: его цель — вычислить одну чисто числовую константу из-за пропорциональности несферичности гравитационного поля выпуклости; масса; Постоянная Ньютона очевидна. Последнее число, которое печатает код, является средним значением$1/r$. Если в начале выбрано {1.1, 0, 0} вместо {0, 0, 1.1}, программа генерирует 0,89 вместо 0,94. Это доказывает, что гравитационный потенциал эллипсоида ведет себя как$-GM/R - hg/2$на расстоянии$R$от происхождения, где$h$— локальная высота поверхности относительно идеализированной сферической поверхности.

В приведенном выше коде я выбрал эллипсоид с радиусами (0,9, 0,9, 1), который представляет собой вытянутый сфероид (длинный, палообразный), в отличие от Земли, которая близка к сплюснутому сфероиду (плоскому, дискообразному). Так что пусть вас не смущают некоторые знаки — они работают нормально.

Бонус от Исаака

Мариано С. указал на следующее решение довольно известного автора:

http://books.google.com/books?id=ySYULc7VEwsC&lpg=PP1&dq=principia%20mathematica&pg=PA424#v=onepage&q&f=false

И) Плоскостность. Здесь мы хотели бы аналитически вычислить решение Любоша Мотля в первом порядке по параметру плоскостности$f$,

$$0~<f~:=~1-\frac{b}{a}~\approx~\frac{a}{b}-1 ~\ll~ 1, \tag{1}$$

где$a$и$b$- экваториальный и полярный радиусы Земли соответственно, и$a>b$. ($\approx$отныне будет означать равенство до членов более высокого порядка в$f$.) Мы предполагаем, что Земля представляет собой массивный сплюснутый эллипсоид

$$\left(\frac{x}{a}\right)^2 + \left(\frac{y}{a}\right)^2 +\left(\frac{z}{b}\right)^2 ~\leq~ 1 \tag{2}$$

с равномерной плотностью$\rho$и объем

$$V~=~\frac{4}{3}\pi a^2 b~\approx~\frac{4}{3}\pi b^3(1+2f).\tag{3}$$

Эксцентриситет

$$e~=~\sqrt{(2-f)f}~\approx~ \sqrt{2f}~\ll ~1.\tag{4}$$

II) квадруполь. Будем считать, что гравитация ньютоновская. Что нам нужно вычислить, так это вклад квадрупольного момента в гравитационный потенциал

$$U_4 ~= ~\frac{x^i Q_{ij}x^j}{r^5},\tag{5}$$

По соображениям симметрии одно из главных направлений квадрупольного момента$Q_{ij}=Q_{ji}$должно быть вдоль полярного$z$-ось, а два других главных направления должны иметь одинаковые собственные значения и лежать в экваториальной плоскости.$xy$самолет. Таким образом$U_4$должен быть в форме

$$U_4 ~=~ \frac{Q_a s^2+Q_b c^2}{r^3}, \tag{6}$$

где$Q_a$и$Q_b$– экваториальное и полярное собственные значения соответственно; где мы ввели сокращенное обозначение$s\equiv\sin(\theta)$и$c\equiv\cos(\theta)$; и где$\theta\in[0,\pi]$является полярным углом. Поскольку квадрупольный момент$Q_{ij}$не может внести вклад в закон Гаусса, мы должны потребовать, чтобы лапласиан$\nabla^2 U_4 = 0$обращается в нуль, что приводит к тому, что полярное собственное значение в два раза меньше экваториального собственного значения, $Q_b=-2Q_a$. Другими словами,

$$U_4 ~\sim~ \frac{3c^2-1}{r^3}~=~\frac{2-3s^2}{r^3}.\tag{7}$$

Поэтому достаточно вычислить гравитационный потенциал в точке$P=(0,0,r)$на полярном$z$-ось, где у нас есть азимутальная$U(1)$симметрия$\varphi\to\varphi+\varphi_{0}$.

III) Кольцо. Вычислим вклад в гравитационную потенциальную энергию кольца, параллельного$xy$плоскость и с центром$(0,0,z)$на$z$-ось. Из-за$U(1)$азимутальной симметрии, мы можем сфокусироваться на точке на кольце с$y=0$и$x\geq 0$, и которые лежат в$xz$самолет. Пусть эта точка имеет двумерные полярные координаты$(R,\theta)$. Другими словами, точка имеет трехмерные координаты.$(x,y,z)=(Rs,0,Rc)$. Из уравнения эллипса

$$\left(\frac{x}{a}\right)^2 +\left(\frac{z}{b}\right)^2 ~=~ 1, \tag{8}$$

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

$$R ~\approx~ b (1 + fs^2) \quad \left\{\begin{array}{rcl}~=~ b &\mathrm{for}& \theta~=~0, \\ ~\approx~ a &\mathrm{for}& \theta~=~\frac{\pi}{2} . \end{array}\right.\tag{9}$$

Далее, расстояние$\ell$с кольца в точку$P$дается косинусным соотношением

$$\ell^2 ~=~ r^2 +R^2 - 2rRc ~\approx~ r^2 +b^2(1 + 2fs^2) - 2rb(1 + fs^2)c$$ $$~=~(r^2 +b^2 -2rbc) +2fbs^2(b-rc).\tag{10}$$

Площадь «внешней» поверхности кольца

$$dA ~=~ 2\pi Rs \cdot R d\theta ~\approx~ 2\pi b^2(1 + 2fs^2) (-dc), \tag{11}$$

с бесконечно малой «радиальной» толщиной

$$dR ~\approx~ (1 + fs^2) db.\tag{12}$$

Таким образом, объем кольца

$$d^2V ~=~ dA \ dR ~\approx~ 2\pi b^2db(1 + 3fs^2) (-dc). \tag{13}$$

(В конце концов мы хотим интегрировать по полярному углу$\theta$от$0$к$\pi$. Это соответствует интегрированию по$c\equiv\cos(\theta)$от$1$к$-1$в отрицательном направлении. Поэтому$dc<0$является отрицательным.)

IV) Потенциал. Вклад кольца в гравитационную потенциальную энергию$d^2 U$в точку$P$является

$$d^2 U ~=~ -G \frac{d^2 M}{\ell} ~=~ - G\rho \frac{d^2 V}{\ell} ~\approx~ G\rho\frac{2\pi b^2db(1 + 3fs^2) dc}{\sqrt{(r^2 +b^2 -2rbc) +2fbs^2(b-rc)}}$$ $$~\approx~ G\rho2\pi b^2db\ dc\left(\frac{1 + 3f(1-c^2)}{\sqrt{r^2 +b^2 -2rbc}}- \frac{1}{2} \frac{2fb(1-c^2)(b-rc)}{(r^2 +b^2 -2rbc)^{\frac{3}{2}}}\right).\tag{14}$$

Интегрирование по косинусу$c$полярного угла$\theta\in[0,\pi]$для формирования тонкой оболочки выходов (с помощью программы MAPLE)

$$dU ~=~ \int_{c=1}^{c=-1} d^2 U ~\approx~ -G\rho\frac{4\pi b^2db}{r}\left( (1+3f) -f( 1+ \frac{2b^2}{5r^2}) - \frac{4fb^2}{15r^2}\right)$$ $$~=~ -G\rho\frac{4\pi b^2db}{r}\left( (1+2f) - \frac{2fb^2}{3r^2}\right).\tag{15}$$

Далее мы интегрируем по$b$чтобы получить полярный потенциал массивного эллипсоида

$$U ~=~ \int^{b} dU ~\approx~ -G\rho\frac{4\pi b^3}{3r}\left( (1+2f) - \frac{2fb^2}{5r^2}\right) ~=~U_1+U_{4b}.\tag{16}$$

V) Монополь и квадруполь. Первый срок$U_1$есть, конечно, именно монопольный потенциал

$$U_1(r)~=~-\frac{GM}{r},\tag{17}$$

и второй срок$U_{4b}$- полярный квадрупольный потенциал

$$U_{4b}~\approx~-\frac{2fb^2}{5r^2}U_1.\tag{18}$$

Таким образом, мы знаем, что гравитационный квадрупольный потенциал в произвольной точке (не обязательно на$z$-ось) есть

$$U_{4}(r)~\approx~-\frac{(2-3s^2)fR^2}{5r^2}U_1(r),\tag{19}$$

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

$$U(r) ~\approx~ U_1(r)+U_4(r) ~\approx~ -\frac{GM}{r}\left(1 - \frac{(2-3s^2)fR^2}{5r^2}\right). \tag{20}$$

VI) Поверхность. С этого момента будем рассматривать только точки$P$с

$$r~=~R~\approx~ b (1 + fs^2)\tag{21}$$

на поверхности эллипсовидной Земли. Тогда монопольный потенциал в$P$становится

$$U_1(R)~ \approx~ -(1- fs^2)\frac{GM}{b}, \tag{22}$$

а квадрупольный потенциал в$P$становится

$$U_4(R) ~\approx~ \frac{f(2-3s^2)}{5}\frac{GM}{b} \tag{23}$$

так что полный гравитационный потенциал в$P$становится

$$U(R) ~\approx~ U_1(R) + U_4(R) ~\approx~ -(1-\frac{2f(1+s^2)}{5})\frac{GM}{b}.\tag{24}$$

VII) Обсуждение. Возьмем за точку отсчета северный полюс, т. е. вычтем гравитационный потенциал$U_b(b)$на Северном полюсе. Тогда гравитационный потенциал в$P$становится

$$U(R)- U_b(b) ~\approx~ \frac{2fs^2}{5} \frac{GM}{b}~\approx~ \frac{2}{5}g hs^2,\tag{25}$$

где

$$h~:~=~a-b~=~fa~ \approx~ fb ~>~0\tag{26}$$

- разница экваториального и полярного радиусов, а

$$g~=~ |\vec{\nabla} U(R)| ~=~ \left. \frac{\partial U(r)}{\partial r}\right|_{r=R}+{\cal O}(f) ~=~ \frac{GM}{R^2}+{\cal O}(f).\tag{27}$$

Теперь мы должны добавить центробежный потенциал

$$V ~=~ -\frac{(\omega Rs)^2}{2}. \tag{28}$$

К заказу в$f\ll 1$что мы работаем, мы видим, что общий потенциал$U(R)+V$постоянна (не зависит от точки поверхности$P$), если

$$\frac{2}{5}gh~=~\frac{(\omega R)^2}{2}. \tag{29}$$

Вывод: находим фактор$\color{Red}{\frac{2}{5}}$в отличие от исходного монопольного аргумента Марка Эйхенлауба.

---

Обновление: через полгода после того, как этот ответ был опубликован на Phys.SE, 18 ноября 2011 г., страница Википедии изменила указанное математическое выражение для параметра плоскостности на$\frac{h}{R}\approx f\approx\frac{5}{4}\frac{\omega^2R}{g}$, и теперь полностью согласен с этим ответом.

Были некоторые сомнения по поводу ответа Любоша (который я принял), так что это просто проверка.

Я скопировал метод, описанный Любошем, и нашел разность потенциалов для эллипсоида с разными эксцентриситетами. Конечно же, для сплюснутого сфероида, если вы сделаете расстояние от центра до экватора долей$e$больше, чем расстояние между центром и полюсом, потенциал составляет примерно долю$e/2$на экваторе меньше.

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

Например, для сплюснутого сфероида

$$\left(\frac{x}{1.01}\right)^2 + \left(\frac{y}{1.01}\right)^2 + z^2 < 1$$

среднее значение$1/r$в$(0,0,1)$составляет около 0,996, а в среднем на$(1.01,0,0)$составляет около 0,991.

Код Python ниже (извините за дилетантство)

import random
import math

points = 10000000
e = .01
rad = 1+e
diam = 2*rad

pot= 0
count = 0
for i in range(1,points):
    x = diam*random.random()-rad
    y = diam*random.random()-rad
    z = diam*random.random()-rad
    r = math.sqrt((x-rad)*(x-rad)+y*y+z*z)
    if x*x/(rad*rad)+y*y/(rad*rad) + z*z < 1:
        pot = pot + 1/r
        count = count + 1

print pot/count

pot2 = 0
count = 0
for j in range(1,points):
    x = diam*random.random()-rad
    y = diam*random.random()-rad
    z = diam*random.random()-rad
    r = math.sqrt(x*x+y*y+(z-1.0)*(z-1.0))
    if x*x/(rad*rad)+y*y/(rad*rad) + z*z < 1:
        pot2 = pot2 + 1/r
        count = count + 1

print pot2/count

Начните с невозмущенного гравитационного потенциала для однородной сферы массы M и радиуса R, внутренней и внешней:

$$\phi^0_\mathrm{in} = {-3M \over 2R} + {M\over 2R^3} (x^2 + y^2 + z^2)$$ $$\phi^0_\mathrm{out} = {- M\over r}$$

Добавьте квадрупольное возмущение, вы получите

$$\phi_\mathrm{in} = \phi^0_\mathrm{in} + {\epsilon M\over R^3} D$$ $$\phi_\mathrm{out} = \phi^0_\mathrm{out} + {M\epsilon R^2\over r^5} D$$

$$D = x^2 + y^2 - 2 z^2$$

Масштабные коэффициенты M и R нужны только для того, чтобы$\epsilon$безразмерный, спад$D\over r^5$просто так, что внешнее решение решает уравнение Лапласа, а согласование решений должно гарантировать, что на любом эллипсоиде около сферы радиуса R два решения равны по порядку$\epsilon$. Причина, по которой это работает, заключается в том, что$\phi^0$решения совпадают как по значению, так и по первой производной при x = R, поэтому они остаются согласованными по значению в ведущем порядке, даже если они отклонены от сферы. Приказ$\epsilon$квадрупольные члены равны на сфере и, следовательно, соответствуют ведущему порядку.

Эллипсоид, который я выберу, решает уравнение:

$$r^2 + \delta D = R^2$$

Z-диаметр увеличивается на долю$\delta$, а диаметр x уменьшился на$\delta/2$. Так что отношение полярного радиуса к экваториальному равно$3\delta/2$. К ведущему заказу

$$r = R + {\delta D \over 2R}$$

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

$$d\phi_\mathrm{in} = {M\over R^3} (rdr) + {\epsilon M\over R^3} dD$$ $$d\phi_\mathrm{out} = {M\over r^3} (rdr) + {MR^2\epsilon \over r^5} dD - {5\epsilon R^2 M\over r^7} (rdr)$$

$$rdr = x dx + y dy + z dz$$ $$dD = 2 x dx + 2ydy - 4z dz$$

На первый заказ в$\epsilon$, только первый член второго уравнения изменяется из-за того, что r не является постоянным на эллипсоиде. Специализируясь на поверхности эллипсоида:

$$d\phi_\mathrm{out}|_\mathrm{ellipsoid} = {M\over R^3} (rdr) + {3\delta \over 2 R^5}(rdr) + {\epsilon M \over R^3} dD - {5\epsilon M \over R^5} (rdr)$$

Приравнивая входную и выходную производные, части, пропорциональные$dD$отменить (как и должно быть — тангенциальные производные равны, потому что две функции равны на эллипсоиде). Остальные тоже должны отменить, так что

$${3\over 2} \delta = 5 \epsilon$$

Таким образом, вы находите связь между$\delta$и$\epsilon$. Решение для$\phi_\mathrm{in}$дает

$$\phi_\mathrm{in} + {3M\over 2R} = {M\over 2R^3}( r^2 + {3\over 5} \delta D )$$

Это означает, глядя на уравнение в скобках, что эквипотенциалы на 60% сжаты, как эллипсоид.

Теперь есть условие, что это уравновешивается вращением, а это означает, что эллипсоид является эквипотенциальным, если добавить центробежный потенциал:

$$- {\omega^2\over 2} (x^2 + y^2) = -{\omega^2 \over 3} (x^2 + y^2 + z^2) -{\omega^2\over 6} (x^2 + y^2 - 2z^2)$$

Чтобы сделать$\delta$эквипотенциальность эллипсоида требует, чтобы$\omega^2\over 6$равно оставшемуся${2\over 5} {M\over 2R^2}$, так что, позвонив$M\over R^2$(ускорение свободного падения) под названием «g», и$\omega^2 R$по наименованию «С» (центробежный)

$$\delta = {5\over 6} {C \over g}$$

Фактическая разница экваториального и полярного диаметров находится путем умножения на 3/2 (см. выше):

$${3\over 2} \delta = {5\over 4} {C\over g}$$

вместо наивной оценки${C\over 2g}$. Таким образом, наивная оценка умножается на два с половиной для вращающейся сферы с одинаковой плотностью.

Неоднородный интерьер: примитивная модель

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

Два конкурирующих эффекта: 1. внутреннее распределение более эллиптическое, чем поверхностное, потому что внутреннее решение чувствует, что вся окружающая эллиптическая Земля деформирует его, а дополнительная плотность деформирует его больше. 2. Эллиптичность внутренности подавляется$1/r^3$спад квадрупольного решения уравнения Лапласа, который$1/r^2$быстрее обычного потенциала. Таким образом, хотя внутренняя часть несколько более деформирована, спад более чем компенсирует эффект внутренней дополнительной плотности, который делает Землю более сферической, хотя и ненамного.

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

Точное решение немного сложное, поэтому я начну с дурацкой модели. Это предполагает, что Земля представляет собой однородный эллипсоид с массой M и параметром эллиптичности$\delta$, плюс точечный источник посередине (или сфера, это не имеет значения), что составляет дополнительную массу внутри, массу M'. Внутренний потенциал задается суперпозицией. С центробежным потенциалом:

$$\phi_{int} = - {M'\over r} - {2M\over 3R} + {M\over 2R^3}(r^2 - {3\over 5} \delta D) + {\omega^2\over 2} r^2 - {\omega^2\over 6} D$$

Он имеет схематическую форму сферического плюс квадруполь (включая центробежную силу внутри F и G).

$$\phi_{int} = F(r) + G(r) D$$

Условие, что$\delta$эллипсоид является эквипотенциалом, находится заменой$r$с$R - {\delta D\over 2R}$внутри F(r) и установив D-часть в ноль:

$${F'(R) \delta \over 2R} = G(r)$$

В этом случае вы получите приведенное ниже уравнение, которое сводится к предыдущему случаю, когда$M'=0$:

$${M'\over M+M'}\delta + {M\over M+M'} (\delta - {3\over 5} \delta) = - {C\over 3 g }$$

где$C=\omega^2 R$центробежная сила, а$g= {M+M'\over R^2}$сила тяжести на поверхности. Следует отметить, что сферическая часть центробежного потенциала${\omega^2\over 2} r^2$всегда вносит подчиненный член, пропорциональный$\omega^2\delta$к уравнению и должны быть опущены. Результат

$${3\over 2} \delta = {1\over 2 (1 - {3\over 5} {M\over M+M'}) } {C\over g}$$

Таким образом, если вы выберете M' равным 0,2 M, вы получите правильный ответ, так что дополнительный экваториальный радиус будет в два раза больше наивной величины${C\over 2g}$.

Это говорит о том, что потенциал на поверхности Земли изменен только по сравнению с оценкой однородного эллипсоида путем добавления сферы с 20% общей массы в центре. Это несколько мало, учитывая, что неоднородная плотность внутри содержит около 25% массы Земли (возмущающая масса вдвое превышает плотность на половине радиуса, то есть около 25% от общего количества). Небольшая разница связана с эллиптичностью ядра.

Неоднородная массовая плотность II: точное решение

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

Вы можете рассматривать его как суперпозицию однородных эллипсоидов среднего радиуса s с параметром эллиптичности$\delta(s)$для$0<s<R$увеличивается по мере продвижения к центру. Каждый из них однороден внутри, с массовой плотностью$|\rho'(s)|$где$\rho(s)$- дополнительная плотность Земли на расстоянии s от центра, так что$\rho(R)=0$. Эти эллипсоиды накладываются на эллипсоид однородной плотности с плотностью$\rho_0$равна поверхностной плотности земной коры:

Я учту$\rho(s)$и$\rho_0$известно, так что я также знаю$|\rho'(s)|$, это (отрицательная) производная по s, которая является плотностью эллипсоида, которую вы добавляете к s, и я также знаю:

$$M(r) = \int_0^r 4\pi \rho(s) s^2 ds$$

Количество$M(s)$является${1\over 4\pi}$раз больше дополнительной массы внутри по сравнению с однородной Землей при плотности земной коры. Обратите внимание, что$M(s)$не зависит от эллипсоидальной формы в ведущем порядке, потому что все вложенные эллипсоиды представляют собой квадрупольные возмущения и, следовательно, содержат тот же объем, что и сферы.

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

$$\phi_\mathrm{int}(r) = F(r) + G(r) D + {\omega^2 \over 6} D$$

Вы знаете, что каждый из этих вложенных эллипсоидов является эквипотенциальным

$$F(s - {\delta(s) \over 2s}) D + G(s) D - {\omega^2\over 6} D$$

так что уравнение, требующее, чтобы это был эквипотенциал при любом s, имеет вид

$${\delta(s) F'(s) \over 2s} - G(s) + {\omega^2\over 6} = 0$$

Чтобы найти форму F и G, вы сначала выражаете внутреннее/внешнее решение для однородного эллипсоида с точки зрения плотности$\rho$и радиус R:

$${\phi_\mathrm{int}\over 4\pi} = - {\rho R^2\over 2} + {\rho\over 6} r^2 + {\rho \delta\over 10} D$$

$${\phi_\mathrm{ext}\over 4\pi} = - {\rho R^3 \over 3 r} + {\rho\delta R^5\over 10 r^5} D$$

Вы можете проверить знак и численное значение коэффициентов, используя правило 3/5 для внутренних эквипотенциальных эллипсоидов, раздельное согласование сферических и D возмущений при r=R и размерный анализ. я положил фактор$4\pi$на дне$\phi$так что правая часть решает уравнение Лапласа в свободной форме с постоянными.

Теперь вы можете наложить все эллипсоиды, установив$\delta$на каждом эллипсоиде быть$\delta(s)$, параметр$\rho$на каждом эллипсоиде быть$|\rho'(s)|$, и$R$быть$s$. Я только собираюсь дать внутреннее решение в r (выполняя интегрирование по частям на сферической части, где, как вы знаете, окажется ответ, и отбрасывая некоторую аддитивную константу C):

$${\phi_\mathrm{int}(r)\over 4\pi} - C = {\rho_0\over 6} r^2 + {\rho_0 \delta(R)\over 10} D - {M(r)\over 4\pi r} + {1\over 10r^5} \int_0^r |\rho'(s)| \delta(s) s^5 ds D + {1\over 10} \int_r^R |\rho'(s)|\delta(s) D$$

Первые два члена являются внутренним решением для постоянной плотности.$\rho_0$. Третий член — это полный сферический вклад, такой же, как и в сферически-симметричном случае. Четвертый член представляет собой наложенный внешний потенциал эллипсоидов внутри r, а последний член представляет собой наложенный внутренний потенциал эллипсоидов вне r.

Отсюда вы можете прочитать сферическую и квадрупольную части:

$$F(r) = {\rho_0\over 6} r^2 + {M(r)\over r}$$ $$G(r) = {\rho_0\delta(R)\over 10} + {1\over 10r^5} \int_0^r |\rho'(s) |\delta(s) s^5 ds + {1\over 10} \int_r^R |\rho'(s)|\delta(s)$$

Так что интегральное уравнение для$\delta(s)$утверждает, что$\delta(r)$форма является эквипотенциальной на любой глубине.

$${F'(r)\delta(r)\over 2r} - G(r) + {\omega^2 \over 6} = 0$$

Это уравнение можно решить численно для любого профиля массы внутри, чтобы найти$\delta(R)$. Это сложно сделать вручную, но вы можете получить качественное представление.

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

Сжатая базовая модель

Точное решение слишком сложно для бумажных и карандашных расчетов, но ищите [здесь]( http://www.google.com/imgres?hl=en&client=ubuntu&hs=dhf&sa=X&channel=fs&tbm=isch&prmd=imvns&tbnid=hjMCgNhAjHnRiM:&imgrefurl= http://www.springerimages.com/Images/Geosciences/1-10.1007_978-90-481-8702-7_100-1&docid=ijMBfCAOC1GhEM&imgurl=http://img.springerimages.com/Images/SpringerBooks/BSE%253D5898/BOK %253D978-90-481-8702-7/PRT%253D5/MediaObjects/WATER_978-90-481-8702-7_5_Part_Fig1-100_HTML.jpg&w=300&h=228&ei=ZccgUJCTK8iH6QHEuoHICQ&zoom=1&iact=hc&vpx=210&vpy=153&dur=4872&hovh=182&hovw=240&tx =134&ty=82&sig=108672344460589538944&page=1&tbnh=129&tbnw=170&start=0&ndsp=8&ved=1t:429,r:1,s:0,i:79&biw=729&bih=483 ) ,вы видите, что разумно смоделировать Землю как две концентрические сферы радиуса$R$и$R_1$с общей массой$M$и$M_1$и$\delta$и$\delta_1$.

я возьму

$$R_1 = {R\over 2}$$

и

$$M_1 = {M\over 4}$$

то есть внутренняя сфера имеет диаметр 3000 км и плотность в два раза больше, что примерно соответствует действительности. Наложение потенциалов и нахождение уравнения для$\delta$s (двухточечное усечение интегрального уравнения), вы найдете

$$-\delta + {3\over 5} {M_0\over M_0 + M_1} \delta + {3\over 5} {M_1\over M_0 + M_1} \delta_1 ({R_1\over R})^2 = {C\over 3g}$$

$${M_0 \over M_0 + M_1} (-\delta_1 + {3\over 5} \delta) + {M_1 \over M_0 + M_1}( -\delta_1 + {3\over 5} \delta_1) = {C\over 3g}$$

Где

$$g = {M_0+ M_1\over R^2}$$ $$C = \omega^2 R$$

- гравитационная сила и центробежная сила на единицу массы, как обычно. Использование параметров и определение$\epsilon = {3\delta\over 2}$и$\epsilon_1={3\delta_1\over 2}$, находят:

$$- 1.04 \epsilon + .06 \epsilon = {C\over g}$$ $$- 1.76 \epsilon_1 + .96 \epsilon = {C\over g}$$

(это точные десятичные дроби, есть знаменатели 100 и 25). Вычитание двух уравнений дает:

$$\epsilon_1 = {\epsilon\over .91}$$

(все еще точные дроби) Что дает уравнение

$$(-1.04 + {.06\over .91} ) \epsilon = {C\over g}$$

Так что фактор впереди$.974$, вместо наивного 2. Это дает экваториальный диаметр 44,3 км, а не 42,73, что достаточно близко, чтобы модель по существу объяснила все, что вы хотели знать.

Значение$\epsilon_1$также интересно, он говорит вам, что ядро ​​Земли на 9% более эксцентрично, чем внешний эллипсоид самой Земли. Учитывая, что точность модели находится на уровне 3%, это должно быть очень точно.

В этом ответе я представлю основу для использования, а затем включу в нее предыдущие ответы. Позвольте мне суммировать ценности, которые мы здесь имеем. Я буду использовать те же обозначения (насколько это возможно), что и все остальные, и Википедию для сплюснутого сфероида, где$a$- большой, экваториальный, радиус.

• Mark1, способ в вопросе,$2 (a-b) = 21.6 km$
• Mark2, метод в ответе, прошлый ответ раз 5/2 для$2 (a-b) = 54.0 km$

Вот мой подход:

Массу Земли можно принять как комбинацию двух фигур: внутренней сферы с радиусом$b$и краевой объем, который представляет собой сплюснутый сфероид за вычетом внутренней сферы. Земля имеет некоторую среднюю плотность, определяемую$\rho = M/V$, но его можно разделить на два разных типа материалов: ядро ​​и кору. Требуемая общая масса будет тогда требовать, чтобы$M = V_{core} \rho_{core} + V_{edge} \rho_{edge}$, в то время как средняя плотность Земли требует$V = V_{core} + V_{edge}$, ограничение$a$и$b$на одну степень свободы. При написании кода мы можем сказать, что$a$всегда подразумевает$b$в соответствии с$b=3V/(4 \pi a^2)$, плотность земной коры также может означать плотность ядра. Затем ограничение эквипотенциальной поверхности фиксирует другую степень свободы, которую можно использовать для итеративного нахождения значения$a$. Иллюстрация:

Потенциал из внутренней сферы прост. Я напишу это для точки на экваторе и на полюсе в сочетании с другими потенциалами.

$$U_a = G \frac{M_{core}}{a} + U_{a,edge} + \frac{1}{2} \omega^2 a^2$$ $$U_b = G \frac{M_{core}}{b} + U_{b,edge}$$

Очевидно, что самое сложное — вычислить потенциал по смехотворно неправильной форме края. Однако перед этим важно подумать о физических последствиях такого взгляда на проблему. Начнем с того, какие плотности относятся к Земле? Здесь указана средняя плотность, а затем плотность в пределах от 0 до 200 км от поверхности.

$$rho_{Earth} = 5520. \frac{kg}{m^3}$$ $$rho_{crust} = 3400. \frac{kg}{m^3}$$

Когда мы на самом деле решим задачу, мы укажем плотность коры, а затем это будет означать плотность ядра сферы. Является ли этот метод точным? Нет. Главное, что он упускает, это то, что центральная сфера не является сферически-симметричным распределением материи. Теоретически плотности будут расслаиваться в соответствии с постоянными потенциальными линиями. Другими словами, если бы у Земли было ядро ​​высокой плотности, оно тоже было бы сплюснутым сфероидом. По этой причине введение внутренней сферы упускает некоторые детали, но эта модель все еще может быть довольно хорошей.

Реализация этого немного сложна, как указывали другие, но сосредоточение вычислений на краевом объеме очень помогает. Я также использовал метод Монте-Карло, но вычел объем ядра. Другими словами, я выстреливал точками случайным образом в крайний объем настолько эффективно, насколько это было разумно. Для этого я использовал взвешенный метод, и, похоже, он сработал. С 5 миллионами итераций я обнаружил случайные вариации в расчете потенциала при$a$иметь стандартное отклонение около$15 m$эквивалента гравитационного потенциала и около$2 m$эквивалента гравитационного потенциала на$b$. Причина более высокого отклонения в$a$связано с тем, что вблизи него масса больше, а выборка была несмещенной в плоскости yz для расчетов$a$потенциал и несмещенный в плоскости xy для расчетов$b$потенциал. В любом случае, когда я повторял$a$значение, я установил допуск для$100 m$, потому что это должно быть значительно выше, чем статистическая вариация, и этого достаточно для расчета выпуклости. Напомним, это мой метод:

Для расчета потенциала при$a$

• Выберите два случайных значения для$y$и$z$значения между$0$и$a$
• Если эти два значения лежат вне эллипса$\frac{y^2}{a^2}+\frac{z^2}{b^2}<1$затем попробуйте два новых значения и повторите попытку (дает около 30% потери эффективности)
• Попробуйте значение x между поверхностью внутренней сферы и внешним сфероидом. Если пара (y,z) не попадает во внутреннюю сферу, выборка между$0$и внешняя поверхность сфероида.
• Подсчитайте потенциал между выбранной точкой (x, y, z) и (a, 0, 0), используя массу M_edge.
• Повторить между (-x,y,z) и (a,0,0)
• Подсчитайте вес этого образца как расстояние между двумя поверхностями, умноженное на два.
• Потенциал в точке (a, 0, 0) представляет собой общий подсчитанный потенциал, деленный на общий подсчитанный вес.
• Повторите аналогичный метод, чтобы найти потенциал в (0,0,b)
• Численно найдите корень, удовлетворяющий условию эквипотенциальности, обсужденному выше.

Я так и сделал, и для разных значений плотности корки получил следующее.

• $\rho_{crust} = 0$, подразумевая, что гравитационное поле нечувствительно к плоскостности, получает$2(a-b)= 21.8 \pm 0.1 km$
• $\rho_{crust} = 3400. \frac{kg}{m^3}$, разумное значение для плотности земной коры, получает$2(a-b) = 34.7 \pm 0.1 km$
• $\rho_{crust} = 5520. \frac{kg}{m^3}$, получается полностью однородная плотность Земли$2(a-b) = 54.0 \pm 0.1 km$

Я подумал, что это хорошие результаты, потому что первый и последний из них близки к предыдущим ответам в пределах числовой погрешности, а разумное значение плотности земной коры приближается к реальному значению$42 km$.

Если кому-то интересно, я могу посмотреть код для этого на github или что-то подобное. В противном случае он немного длиннее, чем другие, размещенные здесь, поэтому пока я не буду загромождать это пространство.

Я смотрю на две модели "жирной земли":

1. сферически-симметричная внутренняя часть с асферическим поверхностным слоем, находящимся в гидростатическом равновесии. Этот анализ обобщает постоянную плотность, принятую в других ответах, и тем самым демонстрирует чувствительность выравнивания к поверхностной плотности. Я сравниваю результат с результатами различных других ответов.
2. Чтобы оценить эффект внутренней асферичности, я снова рассмотрю случай, проанализированный Роном Маймоном, двух концентрических сплюснутых сфер постоянной плотности, находящихся в гидростатическом равновесии. Мой расчет показывает, что влияние сплюснутости ядра на уплощение поверхности мало по сравнению с влиянием плотного ядра на среднюю плотность.

Для обоих анализов я использую следующий результат: рассмотрим тонкую оболочку из материала с плотностью$\rho_s$, средний радиус$r_o$и толщина

$$\delta r_s(\theta)=-\frac{2}{3}f r_o P_2(\cos \theta)$$

(Это полярные координаты с$P_2 = \frac{1}{2} ( 3\cos^2 \theta -1)$, второй полином Лежандра.) Заметим, что$\Delta r_s = \delta r_s(\pi/2) - \delta r_s(0) = f r_0$, так$f$является уплощение. Легко показать, что средний радиус действительно$r_o$а чистая масса раковины равна 0. (Подумайте о поверхностной плотности массы$\rho_s \delta r_s$, что отрицательно, где$\delta r_s<0$. Не обращайте внимания на афизичность отрицательной массы; на практике эта оболочка будет наложена на сферу.)

Тогда в первом порядке по f поле, создаваемое этой оболочкой, равно

$$r > r_o: \phi_e(r,\theta)= \frac{2}{5} f \frac{G \rho_s V(r_o)}{r} \left( \frac{r_o}{r} \right)^2 P_2(\cos\theta)$$

$$r < r_o: \phi_e(r,\theta)= \frac{2}{5} f \frac{G \rho_s V(r_o)}{r_o} \left( \frac{r}{r_o} \right)^2 P_2(\cos\theta)$$

(Здесь$V(r) = \frac{4 \pi}{3} r^3$это объем сферы с радиусом$r$.) Я включаю схему расчета в конце.)

---

1. сферически-симметричный внутренний слой + поверхностный слой

Пусть средняя плотность внутренней части$\rho_a=M/V(r_o)$и плотность на поверхности$\rho_s$. Тогда для описанного выше асферического поверхностного слоя потенциал на поверхности равен (в первом порядке по f):

$$\phi_s = - \frac{GM}{r_{o}}\left( 1+ \frac{2}{3} f P_2(\cos \theta) \right) + \frac{2}{5} f \frac{G V(r_{o}) \rho_s}{r_{o}} P_2(\cos \theta) \\ - \frac{1}{3} \omega^2 r_{o}^2 \left( 1 - P_2(\cos \theta) \right)$$

где первый член включает изменение сферического потенциала («$mgh$"член) второй - квадрупольный потенциал, а последний - псевдопотенциал вращения (записанный в терминах полиномов Лежандра). Грубо говоря, у полюсов поверхность ближе к центру, а значит, и глубже в потенциале хорошо, но уменьшение близлежащей массы частично нейтрализует этот эффект.

Для гидростатического равновесия этот потенциал не должен зависеть от$\theta$, так:

$$-\frac{2}{3} \frac{M}{r_o} f + \frac{2}{5} \frac{\rho_s V(r_o)}{r_o} f = -\frac{1}{3} \omega^2 r_o^2$$

Решение:

$$\left( 1 - \frac{3}{5} \frac{\rho_s}{\rho_a} \right) f = \frac{1}{2} \frac{\omega^2 r_o^2}{GM/r_o} = \frac{1}{2} \frac{\omega^2 r_o}{g} = \frac{q}{2}$$

или

$$\Delta r_s = f r_o = \frac{1}{1-\frac{3}{5} \frac{\rho_s}{\rho_a}} \frac{q}{2} r_o$$

Числа: я использую$r_o=6370$км,$\rho_a = 5520$кг/м$^3$, и$q=3.438 \cdot 10^{-3}$

Случаи.

1. предел 0 поверхностной плотности$(\rho_s/\rho_a \rightarrow 0)$.
   $\Delta r_s = (q/2)r_o = 10.95 \text{ km}$
   Отсутствие поверхностной плотности означает отсутствие квадрупольного поля, поэтому этот результат такой же, как и результат, рассчитанный в вопросе, а также такой же, как в тексте профессора Морина 2004 г., раздел 9.4, задача 8. (Профессор Морин присвоил этой задаче 3 из 4 звезд сложности. ; возможно, это должно было быть 4 из 4?)
2. $\rho_s=\rho_a$, который включает случай постоянной плотности, проанализированный Qmechanic и Роном Маймоном.
   $\Delta r_s = (q/2)r_o /(1-3/5) = (5/2) (q/2) r_o = 27.37 \text{ km}$, в соглашении.
3. $\rho_s=3400$кг/м$^3$, как в численном расчете AlanSE.
   $\Delta r_s = (q/2)r_o /[1-(3/5)(3400/5520)] = 1.586 (q/2) r_o = 17.37 \text{ km}$, опять же согласен.
4. $\rho_s/\rho_a = 4/5$, что соответствует дополнительной массе ядра Рона Маймона.
   $\Delta r_s = (q/2)r_o /[1-(3/5)(4/5)] = (25/13) (q/2) r_o = 1.923 (q/2) r_o =21.06 \text{ km}$.

   ---

2. Две наложенные друг на друга сплюснутые сферы, каждая с постоянной плотностью.

Следуя анализу Рона Маймона, теперь я добавляю основную массу. Обозначение:

1. Я назову сферу меньшего радиуса №1 (как у Рона Маймона) с номинальным радиусом$r_{co}$, масса$M_1$, постоянная плотность$\rho_1=M_1/V(r_{co})$, и сглаживание$f_c$.
2. Сфера большего радиуса - № 2 (отличается от анализа Рона Маймона, извините) с номинальным радиусом.$r_{so}$, масса$M_2$, постоянная плотность$\rho_2=M_2/V(r_{so})$, и сглаживание$f_s$

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

Накладывая поля для двух масс, потенциал на поверхности ядра равен:

$$\phi_c = - \frac{G M_1}{r_{co}} \left( 1+ \frac{2}{3} f_c P_2(\cos \theta) \right) + \frac{2}{5} \frac{G M_1}{r_{co}} \left( \frac{r_{co}}{r_{co}} \right)^2 f_c P_2(\cos \theta) -\frac{3}{2} \frac{G M_2}{r_{so}} \left[ 1 - \frac{1}{3} \left( \frac{r_{co}}{r_{so}} \right)^2 \left( 1 - \frac{4}{3} f_c P_2(\cos \theta) \right) \right] + \frac{2}{5} \frac{G M_2}{r_{so}} \left( \frac{r_{co}}{r_{so}} \right)^2 f_s P_2(\cos \theta) - \frac{1}{3} \omega^2 r_{co}^2 \left( 1 - P_2(\cos \theta) \right)$$

Потенциал на внешней поверхности равен:

$$\phi_s = - \frac{G M_1}{r_{so}} \left( 1+ \frac{2}{3} f_s P_2(\cos \theta) \right) + \frac{2}{5} \frac{G M_1}{r_{so}} \left( \frac{r_{co}}{r_{so}} \right)^2 f_c \, P_2(\cos \theta) - \frac{G M_2}{r_{so}} \left( 1+ \frac{2}{3} f_s P_2(\cos \theta) \right) + \frac{2}{5} \frac{G M_2}{r_{so}} \left( \frac{r_{so}}{r_{so}} \right)^2 f_s \, P_2(\cos \theta) - \frac{1}{3} \omega^2 r_{so}^2 \left( 1 - P_2(\cos \theta) \right)$$

Требование для гидростатического равновесия, чтобы не было$\theta$-зависимость дает пару уравнений:

$$\left[ 1 - \frac{3}{5} \frac{M_2}{M_1+M_2} \right] f_s - \alpha f_c = \frac{q}{2} \text{ , surface}$$

$$-\frac{3}{5} \frac{M_2}{M_1+M_2} f_s + \beta f_c = \frac{q}{2} \text{ , core}$$

$$\text{where } \alpha = \frac{3}{5} \left( \frac{r_{co}}{r_{so}} \right)^2 \frac{M_1}{M_1+M_2}$$

$$\text{and } \beta = \left( \frac{r_{so}}{r_{co}} \right)^3 \left[ \frac{M_1 + \left( \frac{r_{co}}{r_{so}} \right)^3 M_2}{M_1+M_2} - \frac{3}{5} \frac{M_1}{M_1+M_2} \right]$$

Решение:

$$f_c = \frac{1}{\beta + \alpha} f_s$$ $$f_s = \frac{1}{1-\frac{3}{5} \frac{M_2}{M_1+M_2} - \frac{1}{1+\beta/\alpha}} \frac{q}{2}$$

Обратите внимание, что$\frac{M_2}{M_1+M_2} = \frac{\rho_s}{\rho_a}$, в соответствии с предыдущим анализом.

Числа. Используя значения Рона Маймона:$\frac{r_{co}}{r_{so}}=\frac{1}{2}, and \frac{M_1}{M_2}=\frac{1}{4}$, один вычисляет$\alpha = 0.030, \beta = 1.440$и$1/(1+\beta/\alpha) = 0.020$, в результате чего значение для$f_s/(q/2)$рассчитано в части 1 случае 4,$1.923$, увеличивается до 2,002 (+4%), если учитывать сжатие сердцевины. Этот эффект меньше, чем эффект массы ядра в снижении$\rho_s/\rho_a$, что уменьшило$f_s/(q/2)$с 2,5 до 1,923 (-23%).

---

Наконец, вот схема расчета потенциала, связанного с «квадрупольной» оболочкой: аппроксимируйте оболочку как поверхность переменной плотности массы на среднем радиусе$r_o$. Эту поверхность можно разложить на кольца постоянных$\theta$с массой

$$dm(\theta') = 2 \pi r_0 \sin \theta' \, r_0 \, d\theta' \, \delta r_s \, \rho_s = - f \rho_s \frac{4 \pi}{3} r_o^3 \, P_2(\cos \theta') \, \sin \theta' \, d \theta'$$

Потенциал, создаваемый этим кольцом, рассчитан в Jackson, Classical EM , Section 3.3:

$$r>r_o: d\phi (r,\theta) = - \frac{G \, dm(\theta')}{r} \sum_{l=0}^{\infty} \left( \frac{r_o}{r} \right)^l P_l(\cos \theta') P_l(\cos \theta)$$

Подставляя массовую плотность (которая пропорциональна$P_2$) и интегрирование по$\theta'$обнуляет все члены, кроме$P_2$член, в силу ортогональности полиномов Лежандра, дающий указанный результат. (Разработка для$r<r_0$похож.)

Здесь я хотел бы численно проверить теоретическое предсказание фактора$\frac{2}{5}$в отличие от исходного монопольного аргумента Марка Эйхенлауба. На практике это означает вычисление разницы гравитационных потенциалов между Северным полюсом и экватором и деление на соответствующую разницу потенциалов монополей. По численным причинам на практике лучше вычислить обратную (= обратную) дробь, которую затем следует сравнить с$\frac{5}{2}$. Поскольку мои навыки программирования ограничены, я только что написал медленный код MAPLE, чтобы выполнить эту работу.

 b:=100; f:=.10; a:=b*(1+f); V1:=evalf(4*Pi*a^2*b/3); 
 xa:=a; ya:= 0; za:=0; xb:=0; yb:=0; zb:=b; 
 U1a := evalf(V1/sqrt(xa^2 + ya^2 + za^2)); 
 U1b := evalf(V1/sqrt(xb^2 + yb^2 + zb^2));

 Ua:=0;Ub:=0;V:=0;
 for x from -a-.5 by 1 to a+.5 do 
 for y from -a-.5 by 1 to a+.5 do 
 for z from -b-.5 by 1 to b+.5 do 
 if (x/a)^2 + (y/a)^2 + (z/b)^2 < 1 then 
 Ua:=Ua + 1/sqrt((x-xa)^2 + (y-ya)^2 + (z-za)^2); 
 Ub:=Ub + 1/sqrt((x-xb)^2 + (y-yb)^2 + (z-zb)^2); 
 V:=V+1;
 end if;od;od;od;

 b;f;Ua;U1a;Ub;U1b;V;V1;Ub-Ua;U1b-U1a;(U1b-U1a)/(Ub-Ua);

Результат с полярным радиусом$b:=100$был

$$\begin{array}{cc} \mathrm{Flatness}&\mathrm{Inverse \ fraction} \\ f & \frac{U_{1b}-U_{1a}}{U_b-U_a} \\ \\ 0.10 & 2.5549 \\ 0.05 & 2.5353 \\ 0.03 & 2.5149 \\ 0.02 & 2.5085 \\ 0.01 & 2.5088 \\ \end{array}$$

Тот факт, что оценка не улучшается от плоскостности$0.02$к$0.01$является артефактом решетки, потому что шаг решетки того же порядка, что и разница в радиусе между северным полюсом и экватором.