CC BY
40
Вестник Челябинского государственного университета. 2010. № 24 (205). Физика. Вып. 8. С. 64-67.
теплофизика
В. В. Гаврилов
стационарная геотермическая модель
Изучается математическая модель теплового состояния поверхностных слоев Земли в стационарном приближении применительно к району Южного Урала с целью объяснить природу наблюдаемых на поверхности локальных аномалий теплового потока. Восстановление недостающего граничного условия, а также уточнение правой части эллиптического уравнения по данным граничных наблюдений проводится путем решения соответствующих вариационных задач .
Ключевые слова: теплопроводность, температурное поле Земли, тепловой поток, вариационная задача.
Введение. Тепловое состояние Земли характеризуется распределением температуры в ее недрах и теплового потока из недр на поверхности [1] . Плотность теплового потока определяется согласно закону Фурье: q = —X grad T, где Т — температура, X — коэффициент теплопроводности . Тепловой поток через поверхность Земли — важнейшее из всех наблюдаемых геофизических и геологических явлений с энергетической точки зрения, поскольку обусловленная им годовая теплоотдача в десятки раз превышает тепловыделение от сейсмической и вулканической активности [1] .
Хорошо известно существование локальных аномалий теплового потока — зон, в которых его значения резко отличаются по сравнению с сопредельными территориями . Например, на
Урале выделяется весьма протяженная область аномально низкого теплового потока плотностью от 30 мВт/м2 и менее [2] (для континентов эти значения составляют от 20 до 90 мВт/м2) . Проблеме объяснения природы таких аномалий посвящены многие исследования . Применение математического моделирования в геотермических исследованиях ограничено нашими знаниями о глубинном строении Земли и изучаемой территории в частности . Благодаря масштабным исследованиям [2-3] территория Южного Урала в этом отношении изучена достаточно хорошо . Наличие информации о геологическом строении региона (рис . 1) и теплофизических свойствах пород (таблица) позволяет использовать ее при моделировании теплового состояния литосферы
Рис. 1. Модель строения литосферы по профилю «Урсейс»
оценки теплофизических свойств некоторых видов пород
№ Классификация пород Коэффициент теплопроводности, Вт/(м-К) Удельная теплогенерация, мкВт/м3
1 Ультраосновные 2,1—3,1 0,15-0,45
2 Гнейсы, сланцы 2,2-3,3 0,5-2,0
3 Высокобарические 2,3-3,3 0, 3 1 0
4 Восточно-Европейская кора 2,0-3,5 0,2-0,6
5 Мантия 3,0-5,0 0,001-0,02
Хорошо известно, что основным источником тепла в земной коре служит наличие примеси радиоактивных изотопов — урана, тория и калия, концентрация которых убывает с глубиной . Удельную мощность внутренних источников тепла также называют теплогенерацией . Исследование теплового режима приповерхностных слоев Земли представляет интерес с позиции ответа на вопрос, является ли распределение мощности тепловых источников в них основной причиной наблюдаемых на поверхности аномалий теплового потока, или же здесь могут иметь место другие факторы .
Восстановление температурного поля. В силу исключительно высокой тепловой инерции земных недр тепловое состояние литосферы близко к стационарному, поэтому здесь рассматривается модель стационарной теплопроводности . Также предполагается, что преобладает кондуктивный перенос тепла . Зависимость коэффициента теплопроводности от температуры моделируется квадратичной функцией Х(Т) = Х0 - аТ + ЬТ2 либо экспоненциальной Х(Т) = X0в~аТ. Эта зависимость экспериментально изучена для многих пород [2] . Эффект теплового излучения и поглощения среды учитывается внесением в коэффициент теплопроводно-
л гт3 гг!3
сти поправки, равной кизл = с1 =—— 1 ,
3а
где о — постоянная Стефана — Больцмана, а — осредненный по спектру коэффициент поглощения среды [4] . Характерное для пород значение с = 10-10 Вт/(м • К4) [5] . Очевидно, что влияние этой поправки значимо при температурах, близких к 1000 К и выше .
Температурное поле в профиле О = {(х1, х2) | 0 < Хр < /р, р = 1, 2} описывается двумерным уравнением эллиптического типа
2
-Z
Р=1 дхр
к(х, и)
ди
дх
= f (х), х eQ, (1)
в;
где и(х) — температура, к(х, и) — обобщенный коэффициент теплопроводности, /(х) — удельная теплогенерация породы (все величины безразмерные, ось Х1 направлена вниз) . Для корректной постановки задачи необходимы граничные условия В качестве условий на верхней границе Г0 = {(х1, х2) | х1 = 0,0 < х2 < 12} может быть задана температура
либо плотность теплового потока
k^ = 9о(x), x єГ0 dn
(3)
(п — внешняя нормаль), поскольку соответствующие измерения имеются1 . Для задания условий на боковых границах используем следующее соображение . Если бы профиль состоял из слоев, границы которых параллельны оси Х2, а условия на верхней и нижней границах были постоянными, производная равнялась бы нулю (фактически задача свелась бы к одномерной) . С учетом этого расширим исходный профиль (рис . 1), дополнив его справа и слева фиктивными областями слоистой структуры . Тогда на границах расширенного профиля можем задать условия
ди
дП = 0, 0 < < ¡1, х2 = 0, ¡2. (4)
Сложнее обстоит дело с условиями на нижней границе Г, которая недоступна для измерений . Заметим, что на верхней границе Г о мы имеем сразу два условия: (2), (3) . Это обстоятельство позволяет нам сформулировать вариационную задачу для восстановления недостающего граничного условия . Предположим, что мы имеем возможность управлять тепловым потоком внизу:
ди
к — = V(х), х є Г, V є Н, (5)
дп
наблюдая за тепловым потоком на поверхности Функционал качества возьмем в виде
(V) = | С к Щ (Х; V) - #о( Х) ] <ІХ +
+
(6)
а
где + Г, е — параметр штрафа . Второе слагаемое вводится, для того чтобы функционал был сильно выпуклым . Оптимальное управление w (х) определяется из условия
J О) = inf J О) (7)
veH
u(х) = u0(х), x є Г0
(2)
1 Температура на глубине 50-100 м практически по -стоянна и равна примерно 5 °С . Для определения теплового потока существует специальная, довольно сложная методика [1-2], учитывающая многовековое изменение климата, которое существенно искажает измерения на глубинах до 5 км .
или
grad J (w) = 0.
Задача (1), (2), (4)-(7) решается численно на прямоугольной равномерной сетке . Минимизация функционала (6) осуществляется методом градиентного спуска . Для определения неотрицательного оптимального управления используется итерационный процесс:
wn+1 = Wn -Tn+1 Srad J(wn
Wn+1 = max {0, #n+i}, n = 0,1,... (8)
Здесь Tn+1 — итерационный параметр, выбор которого подчинен условию Jg (wn+j) < Je (wn). Если это условие не выполнено, он уменьшается вдвое, а wn+± пересчитывается . Итерации заканчиваются при выполнении условия ||grad Je (w)|| < Z ■ Состояние системы определяется из решения краевой задачи (1), (2), (4), (5) . Для приближенного ее решения строится однородная консервативная разностная схема второго порядка аппроксимации, коэффициенты которой, с учетом разрывности к(x, и), определяются интегро-ин-терполяционным методом [6-7] . Уравнение (1) является квазилинейным . Пусть некоторому изменению управления 8v соответствует изменение состояния бы, причем на интервале [u; и + бы) коэффициенты остаются неизменными для достаточно малого Ъи. Тогда градиент функционала (6) можно выразить через сопряженное
состояние p(х): grad Je (v) = 2| p(x) + \v(x)
V e
Сопряженное состояние определяется на основе разностной формулы Грина, граничных условий (2), (4), (5) и находится из решения краевой задачи, аналогичной (1), (2), (4), (5) . Таким образом, на каждом итерационном шаге (8) решаются две эллиптические краевые задачи по определению основного и сопряженного состояний (используется попеременно-треугольный метод приближенной факторизации [6]) и пересчитываются коэффициенты
На рис . 2 изображено температурное поле, рассчитанное по результатам решения задачи на сетке © = {Ху = (¡Нь у^)} при 0 < і < 54, -15 < j < 200, Ні = 0,115, ^2 = 0,25 (единица соответствует 10 км) . В расчетах использовались осред-ненные теплофизические параметры пород (табл . 1), экспоненциальная зависимость коэффициента теплопроводности от температуры с параметром а = 0,6-0,8 мК-1, в (6) значение параметра е = 200. Точность итерационного приближения принята равной £ = 10-4 . На рис . 3 показано, как выглядит функция теплового потока на поверхности #(х) при соответствующем управлении у(х)|г.
II о I!
Уточнение распределения мощности источников по граничному наблюдению. Исследуем теперь, как влияет распределение мощности тепловых источников в приповерхностных слоях на конфигурацию теплового потока через поверхность . Для этого сформулируем соответствующую вариационную задачу. Пусть состояние системы вместо (1) описывается уравнением
-I
в=1 дхв
к (х, и)
ди
дх.
в
= f (х) + v(х), х eQ,
м
v(х) = Е Vmnm (V
m=1
1, х
< V < V min — ym — Kmax’
Пт ( х ) =
Q
vm’
0, х £ Q
vm
(9)
Здесь ут — постоянные, подлежащие опреде-
м
лению, = У — область управления,
т=1
примыкающая к верхней границе профиля, причем /Х) = /0 при х е О.У. Вместо граничного усло -вия (5) используется условие
к — = q(х), х є Г. дп
(10)
Функционал качества возьмем в виде (6) . Задача (2), (4), (6), (7), (9), (10) по определению
Рис. 2. Модель распределения температуры (в градусах Цельсия) по профилю «Урсейс»
кусочно-постоянного управления решается аналогично предыдущей задаче, с использованием итерационного процесса (8) . Градиент функционала также выражается через сопряженное состояние . Поле теплогенерации, полученное в ре -зультате решения задачи на той же сетке с точностью 10-4, изображено на рис . 4 . На рис . 5 показано, как меняется вид функции теплового по -
тока на поверхности в зависимости от значения теплогенерации в области 1 (см . рис . 1) .
Анализ результатов и выводы. Для оценки погрешности численного решения используется метод, основанный на сгущении сетки [8] . Выполняется серия расчетов при последовательном уменьшении шага сетки вдвое, устанавливаются сходимость и порядок точности
са
о
о
О
-Ч(у1)
•VI
---------д(\2)
----------д(чв)
Км
Рис. 3. Профиль теплового потока на поверхности q(v) и соответствующая функция управления V на нижней границе (VI — оптимальное управление, V2 и Vз — заданные управления)
Рис. 4. Распределение мощности источников (мкВт/м3), оптимизированное по наблюдениям теплового потока на верхней границе
к
о
т
и
р
п
й
о
в
о
-0,3
-0,5
1,0
Км
а
б
Рис. 5. Профиль теплового потока на поверхности при соответствующем значении теплогенерации
(мкВт/м3) в области 1 (см. рис. 1)
Рис. 6. Локальность влияния возмущения теплового потока (мВт/м2) на нижней границе (линией отмечена область заниженных значений)
приближенного решения и производится его уточнение (экстраполяция)
Согласно расчетам, проведенным на основе предлагаемой стационарной геотермической модели с привлечением данных натурного экспери-мента1, значения температуры на глубине 80 км в исследуемом регионе составляют порядка 700800 °С, значения плотности вертикального теплового потока — от 12 до 20 мВт/м2 .
Решение вариационной задачи по реконструкции мощности тепловых источников вблизи поверхности показало наличие области пониженной теплогенерации, соответствующей области пониженного теплового потока на поверхности
Также выявлено, что вид функции q(x)|г весьма чувствителен к значениям теплогенера-ции приповерхностных пород (рис . 5) . Другими словами, тепловой поток на верхней границе гораздо лучше управляется внутренними источниками тепла, чем тепловым потоком на нижней границе (рис З) В частности, при выравнивании теплогенерации в приповерхностных структурах выравнивается и тепловой поток на поверхности, в то время как локальные возмущения теплового потока на нижней границе, напротив, могут не дать заметного эффекта на поверхности при достаточно большой глубине профиля (рис 6)
Относительно влияния теплопроводности пород на конфигурацию поверхностного теплового потока отметим, что значения коэффициента теплопроводности приповерхностных пород несильно различаются и его выравнивание не приводит к существенному изменению вида q ( x) .
Г о
1 Также проводились расчеты, в которых теплофизические параметры задавались как значения нор-
мально распределенных случайных величин .
Анализ результатов расчетов позволяет сделать вывод, что, с точки зрения используемой стационарной модели геотермики, распределение мощности источников тепла в приповерхностных слоях является основным фактором, определяющим конфигурацию теплового пото -ка на поверхности . В частности, падение теплового потока на поверхности является следствием пониженной генерации тепла под поверхностью .
Список литературы
1 . Жарков, В . Н . Внутреннее строение Земли и планет / В . Н . Жарков . М . : Наука, 1983 . 194 с .
2 . Голованова, И . В . Тепловое поле Южного Урала / И . В . Голованова . М . : Наука, 2005. 188 с .
3 . Глубинное строение и геодинамика Южного Урала (проект «Уралсейс») . Тверь : ГЕРС, 2001 . 286 с .
4 . Исаченко, В . П . Теплопередача / В . П . Исаченко, В . А . Осипова, А . С . Сукомел . М . : Энергия, 1975. 488 с .
5 . Zoth, G. Thermal conductivity / G. Zoth, R. Haenel // Handbook of Terrestrial Heat-Flow Density Determination . Dordrecht : Kluver, 1988 . P. 35-41 .
6 . Самарский, А . А. Вычислительная теплопередача / А. А . Самарский, П . Н . Вабищевич . М . : УРСС, 2003 784 с
7. Гаврилов, В . В . Математическая модель температурного режима участка литосферы / В . В . Гаврилов, Н . М . Шерыхалина // Вестн . Уфим . гос . авиац . техн . ун-та . 2007. Т. 9, № 5 (23) . С . 81-86 .
8 Житников, В П Обоснование методов фильтрации результатов численного эксперимента / В П Житников, Н М Шерыхалина // Вестн Уфим . гос . авиац . техн . ун-та . 2007. Т. 9, № 3 . С . 71 .
