Численное решение обратной задачи для нестационарного уравнения конфекции-диффузии Текст научной статьи по специальности «Математика»

УДК 519.63

ЧИСЛЕННОЕ РЕШЕНИЕ РЕТРОСПЕКТИВНОЙ ОБРАТНОЙ ЗАДАЧИ ДЛЯ НЕСТАЦИОНАРНОГО УРАВНЕНИЯ КОНВЕКЦИИ-ДИФФУЗИИ

В, И, Васильев, О, А. Тихонова

Среди некорректных задач для уравнений эволюционного типа особый интерес представляют задачи с несамосопряженными операторами [1]. Примером задач такого типа служит ретроспективная обратная задача конвекции-диффузии.

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

Анализ исследований, в которых рассматривается моделирование процесса конвекции-диффузии, показывает, что не уделяется должного внимания вопросам численного решения ретроспективной обратной задачи конвекции-диффузии. Исключением является работа [1], в которой проводится построение рустойчнвых разностных схем.

Данная работа может рассматриваться как продолжение исследований, начатых в [3,4], где поднимались вопросы численного решения ретроспективной обратной задачи теплопроводности.

1. Постановка задачи

В двумерной области Л, представляющей собой прямоугольник П = {х | х = (х, х), 0 < ха < 1а, а = 1, 2}, © 2009 Васильев В. И., Тихонова О. А.

рассматривается уравнение конвективно-диффузионного переноса:

^ + = хеп, <к*<т, (1)

где слагаемое "У (у)п определяет конвективный перенос.

Для конвективного члена уравнения (1) возможны представления в различных формах: дивергентной, недивергентной и симметричной [2]. Для первого вида записи имеем

= ^^ = vgradw.

ОХ а

а=1

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

\ ^ дуа(х)ш . .

= 2_^ ° = (11У(У«;).

дха

а=1

Обе приведенные формы записи эквивалентны, если движущаяся среда является несжимаемой:

сИУ V = 0.

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

„ 1 / , „ дш дупш\ . .

+ (2) а=1 4 у

Уравнение (1) дополним однородными граничными условиями первого рода:

м(х,г) = о, хе дП, о <г<т. (з)

Требуется найти функцию если начальное условие задано на

момент времени Ь = Т:

и(х,Т) = <р(х), хеП. (4)

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

В связи с этим рассмотрим более подробно оператор конвективного переноса. Для этой цели введем гильбертово пространство Н = и рассмотрим скалярное произведение

{У{ч)и},и1) = J V«; grad ъис1х = ^ J vgradw2<ix = ^ J сИу(«;2у)<±х п п п

— — J ги2 сИуус£х= — J упи)2с1х. п ап

Положим, что нормальная компонента скорости движения среды на границе области обращается в нуль:

уп = чп = 0, х£ дП.

Тогда получим, что

(Г(ч)т,т) = О

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

Г(ч) =

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

диффузии по начальным данным и правой части [2].

В дальнейшем будем рассматривать симметричный вид конвективного слагаемого.

2. Дифференциально-разностная задача

Проведем дискретизацию по пространству. В области Л введем равномерную по каждому направлению сетку & с постоянными шагами

ha, а = 1, 2, и обозначим через ш множество внутренних узлов: ш = {| — ^ X 2_, X 2 ^ ^ ^^ а — ^ а ^^ а ^ ^ а — 1,2,... , 1 ^ ^^ а — 1 а ^ аа — 1, 2}.

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

(у^) = ^2ywh1h2, ||у|| = (у,у)1/2. (5)

Разностный оператор диффузионного переноса на множестве функций у € Н определим выражением

2

Еу= -к^Уха Ха , X € Ш. (6)

а=1

В данном пространстве оператор Л является самосопряженным [5].

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

1 2 1 2

Су = 2+ 2

а=1 а=1

При этом в пространстве Н разностный оператор С будет кососиммет-ричным [2], т. е.

(Су,у)=0, у € Н. (8)

Таким образом, приходим к следующей дифференциально-разностной задаче:

+ Су + Ау = 0, 0 < г < Т, (9)

аЬ

у{х,Т) = ф(х), х € ш, (10)

с положительным самосопряженным оператором Л и кососимметрич-С

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

Предположим, что вместо обратной задачи конвекции-диффузии (9), (10) рассматривается соответствующая прямая задача с начальным условием

у(х, 0) = «(х), х £ &. (11)

Обозначим через уп разностное решение на момент времени = пт, где т > 0 — шаг по времени, причем примем выполненным условие Мт = Т. Поставим в соответствие дифференциально-разностному

уравнению (9) двухслойную явно-неявную разностную схему [2]: уп уп

--— + Суп + Вуп+1 = 0, п = 0,1,..., N-1, (12)

т

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

У0 = V. (13)

Приведенная разностная схема является безусловно р-устойчпвой [2].

3. Вычислительная реализация явно-неявной схемы

Для нахождения решения на новом временном слое из (12), (13) получим сеточную эллиптическую задачу

(1 + т£)уп+1 = ^п, х еП, (14)

с правой частью

= ( 1- тО)уп. (15)

Запишем сеточную эллиптическую задачу (14) в виде операторного уравнения

= ф (16)

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

Л = Е + тВ. (17)

Принимая во внимание свойства самосопряженности и положительности оператора кондуктивного переноса, имеем

А = А*> 0 (18)

при естественном ограничении т > 0.

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

1»к+1 - 1»к

Vfc+l

Awk+1 =Fn, к = 0,1,..., (19)

где vk+i — итерационный параметр.

При применении метода минимальных невязок Vk+i итерационные параметры вычисляются согласно формуле

= rk = Äwk-Fn, к = 0,1,....

(ATk, Ark)

При выполнении двустороннего операторного неравенства [5]

7i > 0, (20)

для числа итераций n, необходимых для достижения относительной точности е, справедлива оценка

. . Ine . .

ü^üi ) = hTp' ( }

где

1-е , 71

р = ТТТ' £= —• 1 + 4 ъ

Принимая во внимание оценки снизу и сверху для оператора Лапласа

к

A к

Ъ = 1 + г-, ^ = 1 + + (22)

Из (22) вытекает следующая асимптотическая зависимость:

С=- = о(-\к\2). (23)

Ъ \т )

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

4. Итерационный метод

Разрешая уравнение (12) относительно у"+1 для заданного приближения у0 па конечный момент времени (п = М), получим

У" = ^ у0, (24)

где Я — оператор перехода с одного временного слоя на другой. Он имеет вид

Я=(Е + тП)— {Е + тС). (25)

Таким образом, приближенному решению обратной задачи соответствует решение следующего сеточного операторного уравнения:

Ли = ф), А = Яи. (26)

В силу несамосопряженности оператора А при решении уравнения (26) используется прием симметризации уравнения, т. е. вместо исходного уравнения рассматривается симметризованное уравнение

Лу = ф, где А = А* А, ф{х) = А*ф. (27)

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

ик+1 ~ ^к +Аик = ф, к = 0,1,... , (28)

где — итерационные параметры.

Для выбора итерационного параметра необходимо ориентироваться на применение итерационных методов вариационного типа. При

использовании итерационного метода минимальных невязок оптимальное значение sk+i берется исходя из условия минимума нормы невязки и вычисляется по формуле

(Ark,rk)

Sk+1 ~~ 7т-

(Air k ,Ark)

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

1. При заданном приближении vk ревизуется прямая задача Avk с использованием разностной схемы nn+1 — nn

---—h Cyn + Луп+1 = 0, n = 0,l,...,N-l, (29)

т

yk = vk, x е u. (30)

Тем самым производится определение сеточной функции yN на конечный момент времени.

yN

решение соответствующей сопряженной задачи Avk = A*Avk: wn — wn+!

--Cwn + Awn = 0, хеш, n = N - 1, ЛГ,... , 0, (31)

т

wN = yN, x е u. (32)

3. Для определения ф = A* ф также решаем сопряженную задачу

- Cz n + kzn = 0, n = N - 1,N,... , 0, (33)

т

zN = ф, x е u. (34)

4. В соответствии с принятой итерационной схемой (28) производится переопределение начального условия:

vk+i = vk - sk+i (ф - Avk). (35)

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

Sfc+1 = , rk = Avk - ф, к = 0,1,.... (36)

(Ark, Ark)

0.25 0.5 0.75 1

Рис. 1. Решение обратной задачи при 6 = 0.

5. Примеры расчетов

Для одномерного случая расчеты проводились при следующих исходных данных: I = 1, Т = 0.3, V = 1 (в безразмерных переменных). По пространственной переменной вводилась равномерная сетка с шагом К = 0.01, то временной — с шагом т = 0.006.

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

(х -0.2)/0.2, 0.2 <х <0.4; (0.5 - х)/0Л, 0.4 < х <0.5. После решения прямой задачи за исходные данные принималось значение решений для временного слоя п = N.

На рис. 1 показаны приближенное и*(х, 0) и точное и(х,0) = х) решения. Несмотря на пониженную гладкость восстанавливаемой функции и(х, 0), разработанный алгоритм находит решение за 8-10 итераций. Дальнейшее продолжение итерационного процесса не приводит к улучшению процесса восстановления. Поэтому итерационное уточнение решения можно прекратить по «слипанию» процесса.

и(х,0) = ^х) =

Рис. 2. Решение обратной задачи при 6 = 0.1.

Для исследования влияния погрешностей на решение обратной задачи конвекции-диффузии входные данные возмущались следующим образом:

ф$(х) = ф{х) + 6а, х £ ш,

где 6 — погрешность, а — случайные величина, равномерно распределенная на отрезке [—1,1].

Итерационный процесс прекращался при выполнении условия

1Ы1 < 6.

На рис. 2 приведены полученные приближенные решения при 6.

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

На рис. 3-5 представлены расчетные результаты для различных коэффициентов диффузии. Расчеты проводились при значениях к = . , . , . .

выше, чем больше коэффициент Пекле. Но следует учитывать, что

Рис. 3. Решение обратной задачи при к = 0.1.

Рис. 4. Решение обратной задачи при к = 0.01.

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

Рис. 5. Решение обратной задачи при к = 0.003.

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

ЛИТЕРАТУРА

1. Самарский А. А., Вабищевич П. Н. Численные методы решения обратных задач математической физики М.: Едиториал УРСС, 2004.

2. Самарский А. А., Вабищевич П. Н. Численные методы решения задач конвекции-диффузии. М.: Эдиториал УРСС, 1999.

3. Вабищевич П. Н., Васильев В. И., Тихонова О. А. Численное решение задачи оптимизации профиля имплантирования в микроэлектронике: Сб. докл. // Третья междунар. науч. конф. «Идентификация динамических систем и обратные задачи» М.; СПб., 1998. С. 107-113.

4. Самарский А. А., Вабищевич П. Н., Васильев В. И. Итерационное решение ретроспективной обратной задачи теплопроводности // Математическое моделирование. 1997. Т. 9, № 5. С. 119-127.

5. Самарский А. А. Теория разностных схем. М.: Наука, 1977.

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

г. Якутск

4 марта 2009 г.