Разностная схема решения уравнения теплопроводности

5.9.1 Разностная схема и число узлов сетки

Уравнение теплопроводности (5-1) решалось с использованием неявной конечно-разностной схемы, близкой к рассмотренной в (Peacemont and Rtachford 1954), но приспособленной к изменяющимся значениям термофизическиx параметров пород и переменным шагам по времени и глубине (Галушкин, 1990; Galushkin, 1997; Makhous et al., 1997):

(5-25)

где и - эффективные теплопроводности соседних элементарных слоёв и , соответственно, на (n+1)-ом шаге времени:

(5-26)

Тридиагональная система уравнений для искомых значений температур в узлах решетки , дополненная граничными условиями на поверxности и в основании области счёта, решалась методом прогонки (см. например, Press et al., 1986; Самарский, Гулин, 1989). В результате получалось распределение температур на следующем, (n+1)-ом, шаге времени. В системе ГАЛО шаг Dz менялся как непрерывная функция глубины от значений 0.5 - 10 м на поверхности осадочной толщи до 1000 -3000 м в основании области счета на глубине 100-200 км. Программа включала процедуру неоднократного переписывания (нового разбиения) области счёта на шаги по глубине всякий раз как их число превосходило 1000. При этом сохранялся линейный рост шагов Dz с глубиной в пределах каждой литологически однородной свиты и непрерывность функции Dz(z) во всей области счёта 0 £ z £. ZM.

5.9.2 Корректность конечно-рвзностной схемы решения

5 .9.2.1 Общая ситуация

Если бы теплофизические параметры были постоянны по глубине и времени и конвективные члены в уравнении (5-1) отсутствовали, то неявная конечно-разностная сxема (5-25, 5-26) уравнения теплопроводности обладала бы безусловной устойчивостью. Однако, мы видели, что теплофизические параметры заметно меняются с глубиной и временем как функции пористости, литологии пород, истории их погружения, кроме того имеют место эффекты выделения и поглощения тепла плавления, консолидации осадков и т.д. Всё это затрудняет явную оценку сходимости разностной схемы. В таких случаях, корректность применяемой разностной схемы проверяется сравнением численных решений с имеющимися аналитическими и полуаналитическими аналогами, а также сопоставлением решений, полученных с использованием различных шагов Dt и Dz.

В системе ГАЛО для оценки точности разностных схем использовались следующие аналитические решения; отложение осадков на фундамент с идентичными теплофизическими свойствами пород осадкочного слоя и фундамента, эрозия однородного полупространства, остывание полупространства с однородной начальной температурой, распределение температур в однородном полупространстве с периодически меняющейся температурой на поверхности (Карслоу и Егер 1964; Тёркот и Шуберт 1985); а также полуаналитические методы оценки искажения теплового потока и распределения температур при отложении осадков на фундамент с теплофизическими характеристиками, отличающимися для пород осадочного слоя (Гольмшток 1979; 1980; Галушкин и Смирнов, 1987). Рассмотрим некоторые из этих решений подробнее.

5.9.2.2 Отложение однородных осадков на однородный фундамент с идентичными свойствами

Поверxностный тепловой поток q(t) и глубинное распределение температур T(z,t), которые устанавливаются при отложении несжимаемых однородных осадков с постоянной скоростью V на фундамент с аналогичными термическими характеристиками, описываются формулами (Карслоу, Егер 1964):

(5-27)

Здесь - функция вероятности, - безразмерное время, k= K/ r×Cp - термическая диффузия, r - плотность, К - теплопроводность, Ср - теплоёмкость осадков и фундамента, t - время, V - скорость осадконакопления, q_¥ есть тепловой поток на большой глубине (невозмущённый тепловой поток на нижней границе), - безразмерная глубина, z - глубина, Z1 = (Z - P)/2, Z2 = (Z + P)/2, - температура на глубине тепловой волны Ökt, соответствующая тепловому потоку q_¥. Процесс осадконакопления, описываемый формулами (5-27), воспроизводился по схеме моделирования бассейнов в соответствии с рис.1-4 и с использованием разностной схемы (5-25, 5-26). Рассматривалось отложение осадков со скоростью V = 1 км /млн. лет в течении 20 млн. лет на фундамент со следующими термическими характеристиками: К = 0,005 кал/см сек°C, Cр=0,25 кал/г°С, k=0,008 см²/сек. В процессе численного моделирования основание области счёта принималось на глубине 200 км и здесь поддерживалась температура 1200°C. При воспроизведении процесса с отложением на поверхности порции осадков мощностью Dz=20 м каждые 20 тысяч лет значения теплового потока на поверхности и распределения температур на глубине, рассчитанные по схеме (5-25, 5-26), отличались от вычисленных по аналитической формуле (5-27) не более чем на 0,1% во всем интервале времени 0 £ t £ 20 млн. лет и глубин 0 £ z £ 200 км.

5.9.2.3 Эрозия однородного полупространства.

Аналитическое выражение для поверхностного теплового потока и распределения температуры с глубиной при эрозии однородного полупространства с постоянной скоростью V имеет вид (Карслоу, Егер 1964):

(5-28)

Здесь безразмерные параметры p, q_¥, T_¥, Z₁, Z₂, и функция Ф(y) те же, что и в формулах (5-27). Сравнение аналитических результатов (5-28) с расчётами в рамках конечно-разностной схемы (5-25, 5-26) проводилось нами для эрозии с постоянной скоростью V = 500 м / млн. лет (t £ 20 млн. лет) и V = 1000 км / млн. лет (t £ 40 млн. лет) и для значений термической диффузии k= 0.008 см²/сек и 0,0095 см²/сек, соответственно. На нижней границе области счёта, на глубине 200 км, поддерживалась температура T=1200°C. Сравнение показало, что вычисления по разностной схеме с шагами эрозии Dz=10 м за 20 и 10 тысяч лет, соответственно, воспроизводит аналитические значения температур и поверхностного теплового потока c точностью не менее 0,1 %. Однако, тот же анализ показывает, что в выборе шага эрозии Dz следует соблюдать осторожность. Это подтверждается данными табл. 6-5, где сравниваются расчёты теплового потока по формулам (5-28) и разностной схеме (5-25, 5-26), выполненные для равномерной эрозии 10 км слоя пород в течение 20 млн. лет при значении фонового теплового потока q_¥ = 0.3 ЕТП с различными шагами эрозии. Сравнение демонстрирует заметное расхождение в результатах численного расчёта q с шагами Dz = 10 и 50 м. Однако, расчёты с шагами Dz = 20 м отличались от значений, полученных с шагом Dz = 10 м, не более чем на 0,003 ЕТП. Следовательно, необходим контроль с вариациями шагов Dz, чтобы гарантировать необходимую точность вычислений.

Табл. 6-5 Сравнение аналитических (5-28) и вычисленных по схеме (5-25), (5-26) тепловых потоков для эрозии однородного полупространства с постоянной скоростью V=500 м/млн. лет

Замечания: q₁ - аналитические значения (10); q₂ - вычисления по разностной схеме (7) с шагом Dz = 10 м; q₃ - то же с шагом эрозии Dz = 50 м. (1 ЕТП = 1×10^-6 кал/см²сек = 41.868 мвт/м2)

5.9.2.4 Остывание однородного полупространства с однородной начальной температурой

Аналитическое решение для распределения температур и теплового потока с глубиной имеет в этом случае вид (Карслоу, Егер 1964):

(5-29)

где Ф(y) определено выше и Ts - температура пород на большой глубине. Контрольные расчёты показали, что даже при относительно грубом временном шаге Dt = 0.7 млн. лет разностная схема (5-25) воспроизводит аналитические значения теплового потока на поверхности с точностью не менее 1% и температур не менее 3 % для всех рассмотренных интервалов глубин 0 £ z £ 200 км и времени 0 £ t £ 100 млн. лет.

5.9.2.5 Отложение осадков на фундамент с отличающимися термофизическими характеристиками.

Корректность разностной схемы (5-25) и выбора шагов Dz, Dt изучалась также на задаче отложения однородных осадков с теплофизическими характеристиками Кs, r_s, Cps (k_s) и генерацией радиоактивного тепла Аs на однородный фундамент с термофизическими характеристиками К_f, r_f, Cp_f и А_f. Полуаналитическое решение этой задачи предполагает отложение осадков с постоянной скоростью (V) и пренебрегает консолидацией осадочных пород. Тепловой поток на поверхности осадков с учётом вклада тепла от распада радиоактивных элементов в осадочных породах получается в этом случае численным решением либо интегрального уравнения с ядром Вольтерра (Гольмшок 1979,1980) либо соответствующего дифференциального уравнения типа (5-1) (Галушкин, Смирнов, 1987). Существенно, что решение задачи является функцией небольшого числа безразмерных параметров:

(5-30)

Для варианта S=1, когда теплофизические свойства осадков не отличались от свойств фундамента, решение описывалось формулами (5-27). Численное воспроизведение этого процесса с помощью разностной схемы (5-25) показало, что всегда возможен выбор шагов Dz, Dt, когда полученные значения поверхностного теплового потока (в том числе и с учётом вклада от тепла распада радиоактивных элементов отличатлись от вычисляемых в рамках полуаналитической модели (Гольмшток 1979, 1980; Галушкин, Смирнов, 1987) не более чем на 1% во всем интервале допустимых значений параметров 0 £ P £ 0,5 и 0,3 £ S £ 1.6.

5.9.2.6 Однородное полупространство с периодически меняющейся температурой поверхности

И, наконец, проверка корректности разностной схемы (5-25) была проведена для аналитического решения задачи с распределением температуры в однородном полупространстве с периодически меняющейся температурой поверхности z=0. Известно, что при поддержании на поверхности периодически меняющейся температуры (T(0,t)=T₀ + DT×cos(w×t)). распределение температур с глубиной даётся формулой (5-19), а тепловой поток имеет вид (Карслоу Егер 1964; Тёркот, Шуберт 1985):

(5-31)

где w=2×p/t - частота временных вариаций поверхностной температуры и c = k/r×C_p - термическая диффузия пород. При численном воспроизведении аналитических решений (5-19, 5-31) встаёт проблема формирования начальныx условий задачи. Если принять, что T(z,0)=То = сonst, а T(0,t)=T₀ + DT×cos(w×t), то значения температур (5-20) и тепловых потоков (5-29) будут достигнуты лишь для времен значительно превосходящих период колебания поверхностной температуры t=2p/w. Чтобы обойти эту трудность, мы формировали начальное условие для численного решения задачи с помощью формулы (5-19), положив в ней t = 0. Расчёты проводились для типичных значений теплофизических параметров литосферы: Ср=0.25 кал/г°C, r = 3.3 г/см³, K= 0.00825 кал/см×°C×сек, т.е. k=0.01 см²/сек и для периода изменения поверxностной температуры t=2 млн. лет. Численные расчёты, проведённые с использованием разностной схемы (5-25, 5-26) с шагом по времени Dt=20 тысяч лет, дали распределение температур с глубиной, отличавшееся от аналитических не более чем на 1% по всему рассмотренному интервалу времени 0 до 10 млн. лет и глубин от 0 до 30 км. Тепловой поток на поверхности оказался более чувствительным к выбору шага по времени и при значении Dt = 20 тысяч лет различие в величинах аналитического (5-31) и численного потоков достигало 6%, но оно оставалось менее 1% для всех Dt £ 10 тысяч лет.

Подводя итоги изложенному в разделе, можно сказать, что величина шагов Dt, Dz, использованных в сравнительных вычислениях, менялась от варианта к варианту, однако, во всех случаях их адекватный подбор позволял обеспечить совпадение вычисленных температур и тепловых потоков с результатами аналитических и полуаналитических оценок с точностью до десятых долей процента. При этом, в ряде случаев к выбору шагов должны предъявляться довольно жесткие требования (например, при анализе режима эрозии; Галушкин, 1988, 1990; Makhous et al., 1997). Сравнение решений показало также, что использование в некоторых системах моделирования (Doligez et al.,1986; Ungerer et al 1987; Welte and Yaclin 1988) расчётных схем со сравнительно грубыми шагами по глубине и времени оправдано лишь для анализа тепловых режимов бассейнов с очень умеренными скоростями отложения осадков.