ФИЗИКА ЗЕМЛИ, 2008, № 1, с. 23-30
УДК 550.311
О СПЕКТРАЛЬНОМ МЕТОДЕ РЕШЕНИЯ УРАВНЕНИЯ СТОКСА
© 2008 г. В. П. Трубицын1, И. Е. Рогожина1' 2, М. К. Кабан1' 2
1Институт физики Земли им. О.Ю. Шмидта РАН, г. Москва, trub@ifz.ru 2GeoForschungsZentrum, Potsdam Поступила в редакцию 08.08.2007 г.
Метод решения уравнения Стокса для сферической модели мантии методом разложения по сферическим гармоникам был разработан в работе [Хейгер, О'Коннелл, 1979]. Однако этот метод применим лишь для случая, если вязкость зависит только от глубины. В этом случае уравнение Стокса сводится к системе независимых уравнений для каждой гармоники. При латеральных вариациях вязкости уравнение Стокса содержит члены в виде произведения гармоник, и все преимущества метода гармоник теряются. В работе [Цханг, Христенсен, 1993] развит метод возмущений для случая, когда члены с произведением латеральных вариаций вязкости на скорость малы. Эти члены сначала вычисляются по предыдущему приближению как функции координат и затем разлагаются в ряд по сферическим функциям. В результате уравнения для гармоник остаются независимыми. Очевидным преимуществом спектрального метода является простота учета эффектов самогравитации и сжимаемости. В этом методе также частично снимаются трудности, связанные с особенностями на полюсах. Этот метод пока не получил практического применения, возможно из-за того приведенные в работе [Цханг, Христенсен, 1993] уравнения содержат опечатки, не разъясненные в литературе. В настоящей работе дается вывод системы уравнений для спектрально-итерационного метода решения уравнения Стокса, и анализируются опечатки в формулах работы [Цханг, Христенсен, 1993] , которые существенно влияют на результаты расчета.
PACS: 91.10.Kg
ВВЕДЕНИЕ
При напряжениях, поддерживаемых в течение более сотен тыс. лет, вещество мантии течет подобно вязкой жидкости и поэтому его движение может быть описано уравнением Стокса. Плотность в мантии меняется по вертикали и по лате-рали. Эти вариации плотности, кроме давления, обусловлены различием химического состава и температуры. В поле тяжести в вязкой мантии более легкие области всплывают, а более тяжелые опускаются. Без тепловой конвекции вещество мантии со временем разделилось бы на слои различной плотности. Однако недра мантии нагреты выше адиабатической температуры. Поэтому случайно поднявшийся со дна мантии элемент вещества (несмотря на адиабатическое охлаждение за счет снижения давления) остается горячее окружающей мантии, и будет постоянно подниматься вверх. Соответственно холодное вещество опускается вниз, и в мантии устанавливается тепловая конвекции, которая перемешивает вещество и тем самым постоянно возобновляет химическую неоднородность мантии.
Эволюция термохимической конвекции в мантии характеризуется изменяющимися во времени пространственными распределениями конвективных течений, химических компонент температуры. Для этого нужно решать систему уравнений переноса массы химических компонент, переноса
импульса (уравнение Стокса) и переноса тепла. Однако задача упрощается, если нужно рассчитать конвективные течения в мгновенной модели современной Земли. Пространственное трехмерное распределение плотности в современной мантии можно найти по данным сейсмической томографии, поскольку вариации скоростей сейсмических волн в основном определяются изменениями температуры, которые непосредственно связаны с изменениями плотности. Таким образом, структуру конвективных течений в современной мантии можно рассчитать, решая лишь одно уравнение Стокса с заданным трехмерным распределением плотности. Зная мантийные течения, можно определить динамический рельеф Земли и крупномасштабные аномалии гравитационного поля. Сравнивая вычисленное и измеренное гравитационное поле Земли, можно судить о достоверности рассчитанной структуры мантийных течений.
Поскольку радиальные и латеральные вариации вязкости в мантии могут достигать трех и более порядков, то для решения уравнения Стокса, вместо классического метода конечных разностей, были разработаны другие методы, в частности метод конечных элементов, позволяющие учитывать большие скачки вязкости. Одним из наиболее известных в настоящее время кодов для расчета конвекции в мантии является СитКом (Калифорнийский Институт Технологии Конвек-
ция Мантия), CitCom (California Institute Technology Convection Mantle) [Zhong et al., 2000]. Тем не менее, численным методам присущ ряд ограничений, которые не позволяют с достаточной точностью рассчитывать динамический геоид. В частности, они не позволяют учесть эффект самогравитации, особенно существенный на больших длинах волн (2-4 сферические гармоники), которые доминируют в наблюденном поле. Для сферических моделей наиболее естественным является спектральный метод разложения искомых и заданных функций в ряды по сферическим гармоникам. В случае, когда вязкость не зависит от угловых переменных, уравнение Стокса сводится к системе независимых линейных уравнений для коэффициентов разложения для каждой гармоники. Это позволяет получать решение, которое практически эквивалентно аналитическому. Однако в общем случае 3-мерного распределения вязкости в уравнение Стокса входит произведение рядов для вязкости и для скорости, и система уравнений по гармоникам не разделяется.
Способ учета нелинейных членов как поправок, получаемых путем сборки и разборки функций по сферическим гармоникам, впервые был предложен в работе [Orsag, 1971]. В работах [Gable et al., 1989; Christensen U., Harder, 1991] итерационный метод учета нелинейных членов, возникающих при латеральных вариациях вязкости, был использован в декартовых моделях. В работе [Цханг, Христенсен, 1993] аналогичным методом была получена система уравнений для коэффициентов разложения по сферическим гармоникам для скоростей течений, вязких напряжений и гравитационного поля с учетом самогравитации. При этом рекомендовалось вводить новую переменную - произведение искомой скорости на вязкость с целью улучшения сходимости итерационного метода. Несмотря на то, что в данной работе декларировались существенные преимущества спектрального метода решения уравнения Стокса для сферичной мантии с латеральными вариациями вязкости, применение опубликованной в работе [Цханг, Христенсен, 1993] системы уравнений для численных решений не получило распространения.
В настоящей работе используется методика разложения по сферическим гармоникам, развитая в [Цханг, Христенсен, 1993], однако, как показано ниже, полученная система уравнений оказывается отличной от указанной работы. Таким образом, опубликованные в [Цханг, Христенсен, 1993] формулы содержат принципиальные погрешности, и их использование привело бы к ошибочным результатам, что, по-видимому, и послужило причиной практического игнорирования метода, первоначально обещавшего большие преимущества.
СВОЙСТВА сферических гармоник и их производных
Поверхностные сферические гармоники Уш(6, ф) являются решениями угловой части уравнения Лапласа в сферических координатах
-(1/sin0)Э/Э0[sin97 ] -- (1/ sin2 0)Э27/Эф2 = l (l +1) Ylm,
(1)
Запишем их в обозначениях, следуя формуле (5.9) работы [Цханг, Христенсен, 1993]
Ylmi(9, ф) = {[(2- 5m0)(2l +1 )Х
1/2 (2)
х (l - m)! ]/[(l + m)! ]}1/2PM(cos0)Cmi,
где Cmi = 5üCos^) + 5asin^).
Присоединенные функции Лежандра Plm(t) = (1 -- t2)m/2dmP(t)/dtm определяются через полиномы Лежандра P¡ (t) = [1/(2 ll!)]dl(t2 - 1)l/dtl.
Учитывая, что J^ cos2 (шф^ф = J^ sin2 (шф^ф =
= п и J2п dф = 2п, легко проверить, что гармоники Ylmi (0, ф) в форме (2) удовлетворяют следующему
J 7
Е 1т1
условию ортонормировки
2п п
УгтРгтЧ' ЯП6= 0 0 (3)
= 4п(5II0(5тт0(5„0 = Ято5и ^тт^И >
Sm0 = 4П.
Формула (5.13) в работе [Цханг, Христенсен, 1993] соответствующая (3), приведена с опечаткой, т.к. в ней sm0 = 4п(2 - 5т0). Однако дополнительный множитель (2 - 5т0) уже был введен в (2) или в формулу (5.9) работы [Цханг, Христенсен, 1993] именно из условия, чтобы гармоники были нормированы только на 4п.
Разложение произвольной функции Е(т, 6, ф) в ряд по полной системе сферических гармоник имеет вид
^ т=11 = 2
^(т, 6, ф) = £ £ £ ^ ,„.,(т) У1 .„,,(6, ф), (4)
I=1т=01=1
где коэффициенты разложения ЕГтч (т) находятся из условия ортонормировки (3)
2 п п
F
lmi
(r) = [ 1/ím0]J dфJF(r, 0, ф) Ylmi sin0d0. (5)
00
В соответственных формулах, приведенных в работе [Цханг, Христенсен, 1993, с. 546], коэффициенты разложения вида (5) имеют предынте-гральный множитель не [1/(4п)], [1/[4п(2 - 5т0)].
Анализ и преобразования рядов по сферическим гармоникам существенно упрощаются, если использовать производные сферических гармоник. Соленоидальное поле Ц = (Цг, Ц, иф) = игег + + иеее + ифеф можно представить в виде суммы двух независимых полей, полоидального и тороидального поля Ц = + Г;. При этом компоненты поля Ц разлагаются в ряд по сферическим гармоникам Уш, а компоненты иб и иф - по производным от сферических гармоник Уш [СИапёгазекЪаг, 1961], которые определяются соотношениями
Ye = Э Y/Эе, Yф = (1/ sin 6)Э Y/ Эф.
(6)
Yее + cYe + Yфф = - LY,
Yфе + cY<p - Y0ф = 0,
YФФе + 2 CYФФ - Y^ = 0
(7)
Геее тлее Т7е 2 т7е ,-» Т7<
+ cY - Y - c Y - 2cY
фф + y^ = _ lyb
где c = ctg е , L = l(l + 1).
Пары производных сферических гармоник удовлетворяют обобщенному условию ортогональности и нормировки, которые можно проверить путем дифференцирования по частям с учетом (7)
2п п 2п п
| М Y^
7фт sin е^е -
I ^ 7фт
7ет sin е ¿е = о,
2п
I м
Y; .m.Yim sinе ¿е +
(8)
2п
+
| ¿ф|
фф
Y/1 sinе ¿е = smlbirb
m lm
ml^ll'^mm' '
с учетом ортонормировки (8) найдем коэффициент /а. Аналогично находится коэффициент /ъ
2п
f a = ( 1/Sml)| ¿ф|[f lYL + f2 Yфm] SinейЮ ,
0
2п
(10)
д = (i/sml)| ¿ф|[ f 1 yL - f 2YU sine ¿e.
В дальнейшем для краткости будем опускать индекс г и в некоторых промежуточных формулах опускать все индексы сферических гармоник I, т и г.
Путем дифференцирования с учетом (1) легко получить соотношения для вторых и третьих производных сферических гармоник
УРАВНЕНИЕ СТОКСА В СФЕРИЧЕСКИХ ПЕРЕМЕННЫХ
Уравнение непрерывности для несжимаемой
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.