Математическая модель температурного режима участка литосферы Текст научной статьи по специальности «Математика»

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ

УДК 519.6(07)

В. В. ГАВРИЛОВ, Н. М. ШЕРЫХАЛИНА

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ТЕМПЕРАТУРНОГО РЕЖИМА УЧАСТКА ЛИТОСФЕРЫ

Рассматривается возможность изучения теплового состояния участка литосферы методом математического моделирования и вычислительного эксперимента. Ставится нелинейная краевая задача стационарной теплопроводности и предлагается метод решения. Для повышения достоверности оценок погрешности приближенного решения и улучшения точности результатов предложено использование новой эффективной методики. Математическое моделирование; вычислительный эксперимент; численная фильтрация; экстраполяция

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

1. ПОСТАНОВКА ЗАДАЧИ

Установившийся процесс переноса тепла теплопроводностью в кусочно-однородной изотропной среде описываем уравнением эллиптического типа [1]

-¿Д(^.“)|т)=/М. (1)

ж = (хі.х-і) Є О.

где «(ж) — температура; к(х, и) — коэффициент теплопроводности; /(ж) — удельная мощность внутренних источников тепла, не зависящая от температуры (радиогенная теплоге-

т

нерация); — расчетная область

/3=1

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

Будем считать, что коэффициент теплопроводности удовлетворяет условиям

0 < «1 ^ к(х,и) ^ к-2, ж £ О

и имеет конечное число точек разрыва первого рода по направлению при фиксированном жз_а (а = 1,2).

Для уравнения (1) ставится краевая задача со смешанными условиями. Граница дО, области П подразделяется на 4 участка, которые условно будем называть верхней, нижней и боковыми границами области. На верхней границе задается условие Дирихле

и(х) = щ(х), х 6 90о, (2)

т. е. известна температура. На боковых границах задается тепловой поток (условие Неймана)

Зи

к(х, и)—- = до а(х), х € <90а, а =1,2. (3) ап

На нижней границе в самом общем случае может быть задано граничное условие третьего рода

Зи

к(х,и)——Ь о(х)и = до з(ж), ж € сЮз, (4) ап

где а(х) > 0 — коэффициент конвективного теплообмена, — производная по внешней нормали к границе .

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

2. МЕТОД РЕШЕНИЯ

Для численного исследования поставленной задачи сначала перейдем от размерных переменных к безразмерным. В качестве опорных значений можем использовать для координат, для темпера-

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

д

'2

у

' дх а=1

- ди дхп

= 08, ж Є П.

(5)

где 08 =

(Хптх)2/

число Остроградского;

граничные условия примут вид

= о,

дп

Є дП0

а = 1,2,

г ди _

к——І- Ві - и = Кг. ап

(6)

(7)

где Кг =

число Кирпичева, В1 =

¿.пип

число Био.

Далее переводим задачу в пространство сеточных функций . Учитывая разрывность коэффициентов, построим однородную консервативную разностную схему на основе ин-тегро-интерполяционного метода [1].

Можем построить сетку ш, равномерную по каждому направлению, либо неравномерную по одному направлению (обусловлено особенностями задачи). Остановимся на случае равномерной сетки ,

і = ОЖ, І = МЪ, ЬаМа = 1а, а =1,2}. Уравнению (1) ставится в соответствие система разностных уравнений

Е

а=1

Мх)у2а) = ф(х), х€ш.

(8)

Коэффициенты нелинейного разностного оператора представляются в виде

X г

аі(жі, х-2, у) = - к(хі — 0,5/іі, х-2 — 0, 5/г-2,

X

- (у(х 1 - /гь х-2 - /г2) + ?у(жі, ж2))) + +к{хі — 0,5/іі,.ї2 + 0, 5/г-2,

- (?у(жі - /іі,ж2 + /г2) + у(х 1,ж-2)))

о,2(жі, ж2, у) = - к(хі — 0,5/іі, ж2 — 0, 5/г-2,

X

- (?у(жі - /гь х,2 - /г-2) + у{х\, ж2))) +

+ /г(жі + 0,5/гі, ж2 — 0, 5/г-2,

X

- (?у(жі + /іі,ж2 - /г-2) + у{х\,ж2)))

(9)

где — безразмерный коэффициент теп-

лопроводности, нелинейно зависящий от температуры.

Граничные условия второго и третьего рода аппроксимируются со вторым порядком на решениях задачи (5)-(7) с использованием интегро-интерполяционного метода.

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

по каждому направлению, при котором ступенчатая аппроксимация будет адекватной. Для решения краевой задачи (5)-(7) в нерегулярной области используется метод фиктивных областей [1,2], обладающий большой технологичностью.

Сеточная задача решается итерационным методом. Наиболее простой итерационный процесс связан с вычислением коэффициентов сеточного оператора по предыдущей итерации. В качестве итерационного метода решения может быть использован нелинейный аналог одного из классических методов [3].

3. МЕТОДИКА ОЦЕНКИ ПОГРЕШНОСТИ И ПОВЫШЕНИЯ ТОЧНОСТИ ЧИСЛЕННОГО РЕШЕНИЯ

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

Используемый здесь способ оценки погрешности основан на идее численной фильтрации приближенных результатов [4-6]. Обычно зависимость вычислительной погрешности от числа узлов сетки (математическая модель погрешности) может быть представлена суммой нескольких слагаемых, например,

гГ1 — г = с\пГр1 + с-2П~р- + ... +

+ сьп~Рь + А (п), (10)

где г — точное значение; — приближенное значение, полученное при числе узловых точек (или числе слагаемых суммы), равном п;

— неизвестные коэффициенты, которые не зависят от ; — произвольные действительные числа .

В могут входить не представленные в сумме (10) слагаемые степенного вида, остаточный член, погрешность округления и многие другие составляющие, обусловленные как численным методом, так и конкретной программной реализацией.

Говоря о числе узлов, будем считать, что на самом деле — число отрезков разбиения или номер последнего узла, при условии, что нумерация начинается с нуля.

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

Задача численной фильтрации заключается в последовательном устранении степенных слагаемых суммы (10) при сохранении значения константы . Первый вариант метода (аналог экстраполяции по Ричардсону), связанный с фильтрацией приближенных результатов, полученных при последовательном увеличении узлов в целое число раз, изложен в [4-6]. Рассмотрим второй вариант (аналог экстраполяции Эйткена), особенно актуальный при решении вычислительных задач в нерегулярных областях.

При решении нелинейного уравнения

методом простых итераций модель погрешности представляется в виде суммы степенных слагаемых

С-2 Г

•2п

' сь 7

Ьп

А (п), (11)

где 7 = ф(х*), |7| < 1, х* — точное решение.

Пусть . Рассмотрим два значения , , вычисленные при числе итера-

ций, равном и соответственно. Применяя метод фильтрации, получим

ф = агп-1 + (1 — а) гп =

= г + ^ + 1 — С17" + ... +

— + 1 - а^) с,-Уп + ...

I3 )

и, из условия ^ + 1 — а = 0, опуская простые выкладки,

у{з) = у__ _ ^ 1

7--? — 1

(12)

Коэффициент 7 3 приближенно можно определить из соотношения

„и-1) Лз-Ц

п — 1

71 —2

уи-Ц _ Лз-Ц

-77 *п-1

т

/уи

— 7з

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

нению), то невозможно говорить о сколько-нибудь содержательном результате фильтрации, т. е. о повышении точности приближенного решения. Проще говоря, фильтрация будет «работать», пока «работает» модель погрешности.

Следующим вопросом является вопрос о надежности получаемых оценок. Самой надежной оценкой погрешности является раз-

ность между приближенным и точным результатом. Точное значение, разумеется, может быть известно только для модельных задач. Согласно [4], выход может быть найден, если вместо точного использовать так называемое эталонное (приближенное, но более точное по сравнению с проверяемым) значение. Это более точное значение может быть всего в 3 раза точнее, чем проверяемое. Вычислить его можно, пользуясь тем же способом, что и проверяемое.

Разность представляет со-

бой оценку погрешности приближенного зна-

Д1)

/2)

яв-

чения гп. Разность =

ляется оценкой погрешности экстраполиро-

(1)

ванного значения или оценкой погрешности оценки погрешности. Отношение Уп =

д|,1] /д

имеет смысл относительной размытости оценки Д„. Если уп -С 1, (т. е. если относительная размытость мала), у нас есть основание доверять полученной оценке. В [4] сформулирован критерий надежности оценки по относительной размытости

1

и„ =

к + 1

где К — коэффициент запаса надежности, наиболее часто используемое на практике значение равно 2.

4. РЕШЕНИЕ МОДЕЛЬНОЙ ЗАДАЧИ

Рассмотрим тестовый пример установившегося теплового режима в твердом теле прямоугольного сечения кусочно-однородной структуры. Расчетной областью будет квадрат ,

^ 0,4, а = 1,2}, состоящий из двух подобластей, в которых теплофизические характеристики материала различны. Температурный режим описывается уравнением (1). Коэффициент теплопроводности зависит от температуры следующим образом:

к (и) = кое'

-Хи

0 < А < 1,

где — базовое (опорное) значение коэффициента теплопроводности. В рассматриваемом случае

{Мх),/(х)} =

2, (^1Л, 0 < X! < 0,2,

'2-'(Ж1+СК4)->+>> ’ ^ ^ Х1 ^

т. е. область состоит из двух слоев. Граничные условия следующие:

«(0,х-2) = 0, 0 < х-2 < 0,4,

«(0,4, х2) = 1п 1,6, 0 < х,2 < 0,4,

Йу

к—— = 0, 0 < х\ < 0,4, Х2 = {0,0,4}.

дх2

Введем в квадратную сетку с шагами /?1 = Л-2 = Л. Уравнение (1) во внутренних узлах аппроксимируется разностным уравнением (8) с коэффициентами (9). Граничные условия первого рода аппроксимируется непосредственно:

у(0,ж2)=0, 0 < х2 < 0,4,

?/(0,4, х2) = 1п 1,6, 0 ^ х,2 ^ 0,4,

а условия второго рода — на решениях задачи

^а2(х 1,х2

1

Ь-2)Ух-2 ~ 2 (аіУЗ:і)хі = О,

0 < х,\ < 0,4, х,2 = 0,

^а2(хі,Х2)ух,

1

2 (аіУхі)хі = О,

0 < х,\ < 0,4, х,2 = 0,4.

Полученная разностная задача на пятиточечном шаблоне запишется в виде системы уравнений с симметричной матрицей:

1-Л'{'у; + -'т.

1 *Л-1

= ¿,і = 0,ДГ. (13)

Здесь N¡1 = 0,4 и обе части (13) умножены на к2. Система (13) решается итерационным методом.

Перейдем к оценке погрешности численного решения. Получим несколько решений на последовательности сеток при удвоении (от 2 до 32) при требуемой точности итерационного приближения к решению є = 10-16. Для оценки погрешности будем использовать значение ^ в одной из наиболее «проблемных» точек — на линии разрыва к(х), напри-мерУ^/2,0.

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

/ ч _ Г 1п(.ті + 1), 0 ^ х\ < 0,2,

4х) - | Ц2Я1 + 0,8), 0,2 < X! < 0,4,

у нас есть возможность сравнить результаты фильтрации с использованием точного значения и фильтрации с использованием эталонного значения. В качестве эталонного значения будем использовать наиболее точное экстраполированное значение.

Сравнение с эталоном у = 0.1823215567935290

1д N

Сравнение с точным значением у =!п 1.2

1д N

Рис. 1. Уменьшение относительной погрешности результатов с последовательным применением численной фильтрации (множество точек 0 соответствует оценке погрешности решения задачи)

На рис. 1 представлены результаты повторного применения численной фильтрации. Полученные результаты позволили установить:

а) приближенное решение сходится к точному со вторым порядком;

б) фильтрация при сравнении с точным и при сравнении с эталонным значением дает практически одинаковые результаты;

в) применение фильтрации позволило получить результаты, точность которых на несколько порядков выше, чем результатов, полученных в результате численного решения задачи.

Так, в результате четырехкратной фильтрации мы получили значение, относительная погрешность которого составляет порядка , что соответствует использованию

численного метода восьмого порядка точности. Относительная погрешность вычисленного при значения составляет более

.

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

Теперь рассмотрим применение фильтрации к приближенным результатам, рассчитанным при увеличении числа итераций (аналог метода Эйткена, рис. 2).

Число итераций

Рис. 2. Уменьшение относительной погрешности результатов с увеличением числа итераций и последовательным применением численной фильтрации

Из полученных результатов непосредственно следует:

а) первичное применение метода позво-

лило повысить точность приближенного результата до порядка , в то время как погрешность приближенного решения составляет порядка ;

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

Невозможность дальнейшего уменьшения погрешности (менее 10-14) связана с накоплением погрешности округления.

ВЫВОДЫ

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

Разработанная и предложенная в [4-6] методика оценки погрешности и повышения точности отличается крайней простотой и основана на анализе приближенных значений, рассчитанных тем или иным численным методом. Ее применение к рассмотренной модельной задаче нелинейной стационарной теплопроводности позволило не только достоверно оценить погрешность приближенного решения, скорость сходимости метода, но и уточнить приближенный результат на многие порядки, что вряд ли представлялось возможным даже при использовании численных методов высокого порядка точности.

Отметим, что применение представленных методов направлено на повышение надежности получаемых оценок, а получение высокоточных результатов способствует этому.

СПИСОК ЛИТЕРАТУРЫ

1. Самарский, А. А. Вычислительная теплопередача / А. А. Самарский, П. Н. Вабищевич. М.: УРСС, 2003.

2. Вабищевич, П. Н. Метод фиктивных областей в задачах математической физики П. Н. Вабищевич. М.: Изд-во МГУ, 1991.

3. Самарский, А. А. Методы решения сеточных уравнений / А. А. Самарский, Е. С. Николаев. М.: Наука, 1978.

4. Житников, В. П. Применение численной фильтрации как средства уточнения результатов вычисления и оценки погрешности /

B. П. Житников [и др.] // Гидродинамика больших скоростей и численное моделирование - 2006 : третья междунар. летн. науч. шк. 2006.

5. Шерыхалина, Н. М. Численная фильтрация данных, искаженных нерегулярной погрешностью / Н. М. Шерыхалина, А. А. Ошмарин // Вестник УГАТУ. 2006. Т. 8, № 1(17).

6. Шерыхалина, Н. М. Методы обработки результатов численного эксперимента для увеличения их точности и надежности / Н. М. Шерыхалина // Вестник УГАТУ. Сер. «Управление, вычислительная техника и информатика». 2007. Т. 9. № 2 (20).

C. 127-137.

ОБ АВТОРАХ

Гаврилов Владимир Владимирович, асп. каф. комп. ма-тем. Магистр техн. и тех-нол. (УГАТУ, 2004). Готовит

á. дис. в обл. вычислит. теплопередачи, оценки погрешно-V сти, оптимизации вычислит.

♦ "л* \1'1"|. 'Г..

Шерыхалина Наталия Михайловна, доц. каф. компьют. мат. Дипл. инж. (УГАТУ, 1993). Канд. физ.-мат. наук (БГУ, 1996). Иссл. в обл. волновых течений жидкости, уединенных волн, методов оценки погрешности численных результатов.