Восстановление граничного режима в задаче моделирования стационарной теплопроводности Текст научной статьи по специальности «Физика»

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

УДК 519.6(07)

В. В. ГАВРИЛОВ

ВОССТАНОВЛЕНИЕ ГРАНИЧНОГО РЕЖИМА В ЗАДАЧЕ МОДЕЛИРОВАНИЯ СТАЦИОНАРНОЙ ТЕПЛОПРОВОДНОСТИ

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

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

д = -кУТ,

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

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

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

Рассмотрим задачу расчета стационарного теплового состояния профиля кусочнооднородной структуры, ограниченного прямо-

угольной областью О = {х | х = (х1, х2), 0< Хр < /р, в = 1, 2}. Стационарное температурное поле в исследуемом профиле описывается уравнением

2 Э

- I —

р=1 Эхр

(,, ч Эи Л к(х, и)----------

Эх

ЭхР

= Д (х), х = (х1, х2) еО, (1)

где и(х) - температура, к(х, и) - коэффициент теплопроводности, Дх) - удельная теплогенера-ция породы. На верхней границе у0 и на боковых границах уь у2 профиля О известна плотность теплового потока:

Эи

к(х, и)— = (х), х еу,, у = 0,1,2. (2)

Эп

2

Условия на нижней границе у = ЭО \ У у у не

У=0

известны, однако на верхней границе задано дополнительное условие:

и(х) = и0(х), хе у0. (3)

Коэффициент теплопроводности ограничен:

0 < к1 < к (х, и) <к 2, х еО, и может изменяться скачкообразно конечное число раз по хр при фиксированном хз-р (в = 1, 2).

Задача (1)-(3) сводится к эволюционной задаче по пространственной переменной х2, которую заменим на ґ (аналог времени), а вместо хі будем писать х. Поскольку на у0 заданы температура и поток, рассмотрим задачу продолжения температурного поля с у0 во всю область О :

Э (, / ч Эи —I к(х, ґ, и)—

Эt I Э/

- Ьи = - Д (х, ґ), 0 < х < I,

Контактная информация: (347) 273-32-00

(5)

где Lu = —— | k(x, t, u) — |/, T = /2 ,

Эx ^ Эx )

с начальными условиями

u(x, 0) = u0(x), 0 < x < /,

—и

-k—(x, 0) = q0(x), 0 < x < /, (6)

Э/

и граничными условиями

- k—(0, /) = q1(t), 0 < / £ Т , (7)

Эx

k — (/, /) = q2(t), 0 < / £ Т . (8)

Эx

Смешанная задача для эволюционного уравнения второго порядка (4)-(8) (ее также называют задачей Коши) относится к некорректно поставленным, поскольку ее решение неустойчиво по начальным условиям (5), (6) [2, 3]. Условная корректность имеет место в классе ограниченных положительных решений, однако численное решение задачи вызывает трудности ввиду ее чрезвычайной чувствительности к погрешности начальных условий. Проблема заключается в построении устойчивого вычислительного алгоритма, позволяющего получить приближенное решение удовлетворительной точности при неточных входных данных.

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

Для численного решения задачи (4)-(8) аппроксимируем ее на сетке юйт =

= ^ = ж, 1 = 0,Ы, ЫН = /\п = пх, п = 0,К, Кх = т} спользуя интегро-интерполяционный метод, построим однородную консервативную разностную схему [4]. Задаче (4)-(8) поставим в соответствие разностную задачу:

- (а(x, /)yx )x - (Ь(x, /)у-)( = ф(x, /),

(x, /) е Юнх,

у(x, 0) = и0(x), (10)

(9)

- -Ь(х, х)у - 1(аух)х =1 ф(х,0) + , (11)

х 2 2 х (11)

0 < х < I,

-\ф, ,)ух -\(Ьу,), = 2Ф(0,0 + , (12)

0 < , < Т,

Т7а(1’ґ )л- 2(Ьу'х= 2 ^ґ)+«Г • (13)

0 < , < Т,

у(х + И, ґ) - у(х,,) у(х,,) - у(х - И,,)

где ух = —---—-, Ух = ——--------------------- .

ИИ

Коэффициенты а(х, ґ), Ь(х, ґ) определяются ин-

тегро-интерполяционным методом [4, 5]. Таким

образом, мы пришли к трехслойной явной разностной схеме вида

Ьп+1 Уп+1 - (Ьп +1 + Ьп )Уп + ЬпУп-1

-ЛУп =-Фп ,

(14)

п = 1, К -1,

где оператор Л определяется с использованием (9), (12), (13):

Лу =

2---аух,

И

-(аух )х 2

~аух,

И

х = 0,

И < х < I - И, (15)

х = I,

где а = а(x + Н, /).

Запишем (14) в виде двухслойной схемы

У п +1 = Эу п + *„ , 06)

выбрав у п ={Уп-1, Уп }. Получим:

8 =

( 0 —^Е

V Ьп+1

( 0 >

2 Фп

Е

1+

Е+

п +1 )

п +1

Л

(17)

-X

V

Ь

п+1 )

Нетрудно видеть, что |8|| > 1, поэтому погрешность в ходе вычислений по схеме (14) будет накапливаться (схема неустойчива по начальным условиям).

Построим на основе (14) устойчивую схему, пользуясь принципом регуляризации. Запишем трехслойную схему (14) в каноническом виде:

В Уп+1 Уп—1 + х2СУп+1 2Уп + Уп-1 +

2х х2

+ АУп =-Фп , п = I2,-,

(18)

где

А = -Л, В = Ьп+1 Ъп Е, С = Ьп+1 +Ъп Е .

х 2х2

Регуляризацию проведем за счет возмущения оператора С:

С = Ьп+1 + Ьп е + аЯ,

2х2

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

Ьп+1Уп+1 -(Ьп+1 + Ьп )Уп + ЬпУп-1 -

х2

- ЛУп + «^Уп +1 - 2Уп + Уп-1 ) = -Фп, (19)

п = 1, К -1.

2

х

2

Ь

х

п

Выбор

Я = Л*Л = Л2 (20)

соответствует методу квазиобращения [3], широко используемому при решении некорректных задач. Оператор (20) определяется согласно (15) и, таким образом, схема (19), (20) является пятислойной. Вычислительная реализация схемы связана с использованием метода предиктор-корректор. Следует отметить, что затраты на реализацию данной схемы довольно велики, учитывая нелинейность задачи и присутствие параметра регуляризации, для определения значения которого требуется дополнительный итерационный процесс. Это является основным недостатком схемы.

Варианту упрощенной регуляризации соответствует выбор

Я = Л . (21)

В этом случае мы приходим к неявной трехслойной схеме, которая является аналогом схемы с весами. Счет по схеме (19), (21) устойчив при а > 0. Недостатком схемы является отсутствие регуляризатора высокого порядка, способного производить эффект сглаживания.

Рассмотрим возможность применения послойного сглаживания решения с целью уменьшения влияния погрешности на каждом слое. Используем идеею метода предиктор-корректор: сначала получаем решение 'уп+1 на новом слое, а затем проводим его сглаживание с помощью оператора S:

Уп+1 = 5~п+^ п = 0,1,.... (22)

Алгоритм сглаживания сеточной функции v(x) построим на основе вариационной регуляризации А. Н. Тихонова [7]. Некорректную задачу аппроксимации

г(x) = у(x), у(x), г(x) еюН, 0 £ x £ / (23)

заменим задачей минимизации функционала:

\\г - УІ + аТ( г) ® тіп,

(24)

где Т(г) = НіРу(5)

0 I У=0

- тихоновский

стабилизатор порядка т (р() - весовые коэффициенты), а - параметр регуляризации. Уравнение Эйлера для задачи (24) имеет вид:

т Иг

г(5) + а I (-1)г —

г=0 0 < 5 < І,

(

. ,ИГг Рг (5)^Г

= У(5),

(25)

И

Рг (5)

Игг

= 0,

(26)

г = 1, т, s = 0, /.

Остановимся на использовании стабилизатора порядка т = 2 при р2(^) = 1, р0^) = р1(^) = 0. В этом случае сглаженная функция определяется из решения краевой задачи

(аг1¥ + г = у(x), 0 < x < /

{ г*( x) = г *(x) = 0, x = 0, /

(27)

методом прогонки.

Отметим, что вычислительная реализация построенного алгоритма (23), (27) весьма проста, однако проблема выбора параметра регуляризации остается.

Рассмотрим еще один алгоритм сглаживания, без использования параметра регуляризации. Идея заключается в применении сглаживающего (усредняющего) оператора S, построенного на основе квадратурных формул. Используя формулу трапеций для отрезков [xi-1, xi ] и [xi, xi+1 ], определим значение сглаженной функции в узле xi следующим образом:

1 xi+1

Sv(xi) = — I* у(г. )ix =

2Н/

xi-1

= 1 (v(xг-1) + ^) + v(xг+l )), (28)

1 = 1, N -1.

Использование формулы Симпсона для частичного отрезка [xi-1, xi+1 ] дает:

1 x^+1

Sv(xi) = — I* v(x. =

2Н }

xi-1

=1 ) + 4v(xг) + v(xг +1)),

6

(29)

і = 1, N -1.

а соответствующие граничные условия

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

(7), (8).

3. МОДЕЛЬНАЯ ЗАДАЧА

Рассмотрим модельную задачу расчета стационарного теплового состояния кусочнооднородного образца прямоугольного (квадратного) сечения О = {х | х = (х1, х2), 0 < ху < 0,4,

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

т

]~г

г

к (и) = к0в Хи, 0 <1< 1,

Причем

{ко(х), /(х)} =

3,0625(1 + 1)

1,96, (1,25х1 + 1,15)2+х :

0 < х1 < 0,2, 4,2875(1 + 1)

1,4, (1,75 х + 1,05)2+х

0,2 < х < 0,4 .

Граничные условия имеют вид:

, ди 196 „ \

к =-----------+ 6-с(х1),

дх1 115 Ы 17

х1 = 0, 0 < х1 < 0,4,

(30)

к— = 0, 0 < х1 < 0,4, х2 = 0,0,4, (31)

дх2

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

На участке границы {х1 = 0,4, 0 < х2 < 0,4} условия неизвестны, однако имеется дополнительное граничное условие на участке {х1 = 0, 0 < х2 < 0,4}:

и(0, х2 ) = 1п(1,15). (32)

Для оценки погрешности полученного решения применим метод численной фильтрации приближенных результатов [8-10]. Получим несколько приближенных решений на последовательности сеток при удвоении К. В качестве контрольной точки выберем у(0,4, 0,2). По-

скольку для рассматриваемого примера известно точное решение:

у(0,4, 0,2) = 1п(1,75),

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

Графики, представленные на рис. 1-4, показывают, как изменяется относительная погрешность численного решения с последовательным применением метода фильтрации. Множество точек 0 соответствует оценке погрешности численного решения, рассчитанного по разностной схеме, 1 - оценке погрешности экстраполированного решения и т. д. На диаграммах с литерой А показана сходимость приближенного решения к эталонному, с литерой Б - к точному решению. В расчетах использовалось значение

5 = 10-6.

4. ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Как видно из диаграмм, все четыре алгоритма позволили получить решение, имеющее первый порядок точности (об этом свидетельствует наклон линии 0). Численная фильтрация приближенного решения в первых трех случаях позволила повысить порядок точности до второго (рис.1-3, линия 1), а относительную погрешность уменьшить до 10-6 , т. е. до уровня погрешности начального условия (30).

... о

с г'

: .4

V '*-0 Г'

• *

0,5

1,5

1

■€>■■ 2

2,5

С

У . 4 г''

■■■ й .<

• *

0,5

1д К

1 1,5

1д К

0

1

2

2,5

0

0

2

0

2

а

б

Рис. 1. Оценка погрешности решения, полученного по схеме (19), (20)

с 5*

,4

У

0

1 2

0 0,5 1 1,5 2 2,5

1д К

•

...( У

0

1

■-€>-- 2

0 0,5 1 1,5 2 2,5

1д К

б

Рис. 2. Оценка погрешности решения, полученного по схеме (14) с применением алгоритма сглаживания

(23), (27)

,.-0

.-О

с •Г

/ У''

ж б

»- 0

♦ - 1

о-- 2

0 0,5 1 1,5 2 2,5

1д К

..О''

С У

/ . -< У

в 1 ■ ”

•

»- 0

1

о- 2

0 0,5 1 1,5 2 2,5

1д К

б

Рис. 3. Оценка погрешности решения, полученного по схеме (19), (21) с применением

алгоритма сглаживания (22), (28)

а

а

4 / ► /

/ С Г

♦’ 0* .Л У

1 1,5

1д К

0

1

2

/ Р ► /

/ С

' О' . ■ 1 У

0,5

1 1,5

1д К

0

1

■-€>-- 2

2,5

0

2

0

0,5

2

2,5

б

а

Рис. 4. Оценка погрешности решения, полученного по схеме (14) с применением алгоритма сглаживания

(22), (28)

Повторная фильтрация не привела к улучшению точности. Расчет по явной схеме (14) с применением алгоритма послойного вариационного сглаживания (23), (27) позволил получить решение не меньшей точности, чем при использовании более затратной схемы (19), (20) (рис.2). Различия в диаграммах А и Б на рис. 1, 2 объясняются тем, что точность полученного решения не может выше точности начальных условий. Использование схемы с весами (19), (21) совместно с алгоритмом сглаживания (22), (28) также дает хороший результат. Наконец, использование явной схемы (14) совместно с алгоритмом сглаживания (22), (28) является наиболее экономичным вариантом расчета. В этом случае мы также получаем сходящееся решение, однако экстраполировать его не удается (рис. 4).

ВЫВОДЫ

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

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

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

1. Голованова И. В. Тепловое поле Южного Урала. М.: Наука, 2005. 189 с.

2. Адамар Ж. Задача Коши для линейных уравнений с частными производными гиперболического типа. М.: Наука, 1978. 160 с.

3. Латтес Р., Лионс Ж.-Л. Метод квазиобращения и его приложения. М.: Мир, 1970. 280 с.

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

5. Гаврилов В. В., Шерыхалина Н. М. Математическая модель температурного режима участка литосферы // Вестник УГАТУ, 2007. Т.9, №5 (23). С. 81-86.

6. Самарский А. А., Вабищевич П. Н. Численные методы решения обратных задач математической физики. М.: Изд-во ЛКИ, 2007. 480 с.

7. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач. М.: Наука, 1979. 285 с.

8. Житников В. П. и др. Применение численной фильтрации как средства уточнения результатов вычисления и оценки погрешности // Третья междуна-родн. летн. науч. шк. «Гидродинамика больших скоростей и численное моделирование-2006», 2006.

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

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

ОБ АВТОРЕ

Гаврилов Владимир Владимирович, Дипл. магистр техники и технологий (УГАТУ, 2004). Готовит дис. в обл. мат. моделирования, оптимизации вычислительных методов, комплексов программ.