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

Вычислительные технологии

Том 7, № 6, 2002

РАЗНОСТНАЯ СХЕМА ДЛЯ РЕШЕНИЯ КОНВЕКТИВНО-ДИФФУЗИОННЫХ ЗАДАЧ ТЕПЛОМАССООБМЕНА

В. Г. Зверев, В. Д. Гольдин Томский госуниверситет, НИИ прикладной математики и механики при ТГУ, Россия e-mail: [email protected]

The three-point finite-difference monotonous scheme on non-uniform grid for solving convective-diffusion transport equation with source, discontinuous coefficients and boundary conditions of a general type is suggested. The construction of the scheme is based on the application of analytical solution of non-homogeneous differential equation of the second order with constant coefficients and linear right-hand side. This solution is locally exact within the grid step interval. The asymptotic of scheme coefficients for small values of grid parameter is investigated and the relation to spline approximation schemes is shown. The efficiency of difference scheme application is verified by test calculations.

Введение

Конвективно-диффузионный перенос играет определяющую роль во многих процессах тепломассообмена [1]. В качестве базовых математических моделей при его описании выступают краевые задачи "конвекции — диффузии — реакции", отличительным признаком которых является наличие в определяющих уравнениях слагаемых с первой производной и младших членов. Получение решения, даже в одномерном случае, требует привлечения вычислительных средств, поэтому для аппроксимации краевых задач, обладающих достаточной гладкостью, разработаны эффективные численные методы [2-7].

Доминирование конвекции над диффузией, знакопеременность коэффициента при конвективном слагаемом делают возможным появление в расчетной области локальных областей с большими градиентами искомой функции — пограничных и внутренних переходных слоев. Наличие источников и стоков в химически реагирующих потоках вносит дополнительные сложности в картину процесса. В конечном счете, это приводит к серьезным трудностям при численном интегрировании математических моделей [4, 8, 9].

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

© В. Г. Зверев, В. Д. Гольдин, 2002.

разностных схем, в частности, основные способы задания регуляризаторов обсуждаются в [4, 5].

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

Практические потребности решения таких задач дали импульс к их исследованию [10] и разработке специальных вычислительных технологий, обеспечивающих необходимый порядок точности приближенного решения, равномерно по малому параметру [11, 12]. Данной проблеме посвящено большое количество работ отечественных и зарубежных исследователей, обширная библиография по этому вопросу приводится в монографиях [13-16], статьях [17-23].

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

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

1. Математическая постановка задачи и ее анализ

Рассмотрим одномерное стационарное конвективно-диффузионное уравнение переноса общего вида в несамосопряженной форме на отрезке [а, Ь]

Ьи = ¿х ¿с) + ¿ж + /гз(х)и = Х € (а,Ь) (1)

с краевыми условиями

¿и

х = а : Ьги = дг — + <22и = дз, (2)

¿и

х = Ь : Ь2и = 'Рг— + Р2и = рз. (3)

ах

Здесь и(х) — искомая функция аргумента х, коэффициенты уравнения /г(х) — /4(х) принадлежат классу ^0[а, Ь] кусочно-непрерывных на [а, Ь] функций; дг, рг (г = 1 — 3) — заданные числа, причем д\ + > 0, р\ + '2 > 0 и рг, р2, дг > 0, д2 < 0. Считаем, что выполняются обычные для задач тепломассообмена условия /г(х) > е > 0 и /3(х) < 0, значение е может быть сколь угодно малым.

К необходимости решения уравнения (1)-(3) приводят задачи конвективного тепло-массопереноса в химически реагирующих средах [1]. Экспоненциальный характер решения и наличие узких областей с большими градиентами при значениях |/2(х)| ^ Л(х) характерны для этого уравнения.

Данную ситуацию хорошо иллюстрирует простейший пример стационарной конвекции и диффузии [8, 9, 11], дифференциальное уравнение которого имеет вид

^ и _ , _ ч ^ , \

£5? - И = ° х е ^ (4)

х = 0: и = и0, £ = 1/Ре, Ре = ^Ь/Г, х = 1 : и = и1,

где V — скорость; Ь — масштаб длины; Г — коэффициент диффузии; Ре — число Пекле; х — безразмерная координата.

Точное решение и(х) и его производная и'(х) зависят от числа Ре:

и(х) - и^ ехр(Ре х) - 1

-= £(Ре,х) =-——-—, (5)

и1 — и0 ехр(Ре) — 1

и'(х) _ . Ре ехр(Ре х) , ,

п(Ре, х) = --/. (6)

u1 — u0 ' exp(Pe) — 1

Из (5), (6) следует, что при сильной диффузии Pe ^ 0 (е ^ то), решение u(x) стремится к линейной функции, так как

lim £(Pe,x) = x, lim n(Pe, x) = 1.

Pe—ü Pe—ü

При сильной конвекции Pe ^ то (е ^ 0) функция £(Pe,x) близка нулю везде, кроме правой границы. Здесь в узком слое толщиной около е возникает деформация решения с неограниченным ростом производной

lim £(Pe, x) = exp (—Pe(1 — x)) , lim n(Pe,x) = Pe exp (—Pe(1 — x)).

Pe—Pe—

Зависимости u(x) и u'(x) для различных значений Pe показаны на рис. 1. Отрезок [0,1] без потери общности можно считать за шаг грубой сетки. Линейный профиль (Pe = 0), обычно полагаемый между узлами при построении конечно-разностной схемы, не соответствует действительности при малых значениях е (при Pe = 20, кривые 5) и оказывается пригодным лишь до Pe < 2. Поэтому в полуцелом узле возникает значительная ошибка в величине диффузионного потока q = Гм'(0.5), определяющем точность разностной схемы.

На рис. 2 представлены результаты расчета задачи при Pe = 20 по классическим схемам, построенным на предположении о линейности профиля на сеточном шаге. На грубой сетке с шагом h = 0.2 (соответствует сеточному Pe^ = vh/Г = 4) схема с центральными разностями дает значительные колебания, односторонние разности сильно размазывают решение. В этом случае с погрешностью O(h2) аппроксимируется уравнение

/-, „ ^ ,d2u du e(1 + 0.5Peh)-^ — — = 0, dx2 dx

в котором при значении Pe^ > 2 схемная вязкость превышает физическую [9]. Еще более сильное отличие от точного решения наблюдается в окрестности пограничного слоя для градиента функции (рис. 2, б).

О 0.2 0.4 0.6 0.8 Ж 1.0 0 0.2 0.4 0.6 0.8 X 1.0

а б

Рис. 1. Профили функции и(х) (а) и ее производной и'(х) (б) при Ре = 0, 1, 2, 10, 20 — кривые 1 — 5 соответственно.

Рис. 2. Численное решение конвективно-диффузионного уравнения при числе Ре = 20 на основе различных схем: 1 — точное решение, 2 — центральные разности, 3 — противопотоковые разности, кривая со значком • — настоящая работа, к = 0.2, Ре^ = 4.

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

2. Методика построения разностной схемы

Введем на [а, Ь] неравномерную сетку

^ = {ж^, 2 = 0, N Х0 < . . . < < Xi < Хг+1 < . . . < Хм = Ь} ,

Л — х^1 Жi, Л — Жi Жi_1, Л — Xi_1).

i

Будем считать, что точки разрыва функций /^£) — /4(£) совпадают с узлами сетки. Выделим на П трехточечный шаблон {£¿-1, £¿+1}. Обозначим через

Х1+1 Xi

/+ = У /1(*Э, /- = р / (7)

хг хг—1

квадратурные аппроксимации интегралов с погрешностью 0(Л,2). Аналогичные формулы пусть имеют место для функций /2(ж), /з(х). Здесь и в дальнейшем величинам, относящимся к [£¿-1, £¿1 и [£г, £¿+1], будем приписывать верхние индексы "—" и " +" соответственно. Наряду с (1), на ^, £¿+1] рассмотрим близкую задачу с кусочно-постоянными коэффициентами и правой частью /4(£), линейно зависящей от аргумента £:

= "Т (Л+1Г1 + /2+ТГ + Л+и = (£), £ е ^^ а£ \ а£ / а£

¿ • ГГ ,

£ = ^¿+1 : Г = г(£^+1) =

/+ (£) = /+ £) + /+¿+1. (8)

На отрезке [£¿-1, £¿] задача (8) для оператора выглядит аналогичным образом. Общее решение неоднородного дифференцального уравнения второго порядка (8) может быть записано в виде [32]

( ) = /+ [Л+Л+^+1 — £) — (А+ + Л+)] + Г(£) /1+ Ь+ (Л+Л+)2 +

+ % [Л+ Л+(£ — £^ + (Л+ + Л+)] + С1вл+х + ^ (9)

/1+ й+ (Л+Л+)2

где С1, С2 — константы; Л+ и Л+ — корни характеристического уравнения

Л2 + Л + ^3+ = 0,

/1 /1

Л+ = Л+ = В = (/2+)2 — 4/+/+. (10)

2/+ 2 2/:

1

Так как /3(£) < 0, то В > 0, при этом корни имеют разные знаки. При /3(£) = 0 один из корней равен нулю. Если /3(£) = 0, /2(£) = 0, то оба корня равны нулю.

Постоянные С1 и С2 определяются из краевых условий (8). В конечном итоге точное решение уравнения (8) на отрезке £¿+1] запишется в виде

г(£) = г__I_еЛ+^+

г(£) г (л+ — л+) е

+ 'Г + +

-Л+(жг+1-ж) _ е-Л+(жг+1-ж)

е

+ + /+

+ ^+1С+

N + + (М + — N+) еЛ1 с — (М + — N+) еЛ1(х-хг) + [М +/^+(£¿+1 — £) — N+]} +

+

/+

4,г+1

тг

— (Ы + + N+) с — N+ел+(х-х4)+ [ы+/н+(х — хг) + N +] } , (11)

где

+

(Л+ — А+) н+

e-л+h+

Опираясь на выражение (11), вычислим диффузионный поток д+ (без учета знака) на [хг, в точке х = хг+0 (справа):

Я

+ /+

1 ^х

«¿+0

игЛ+( Л+ — + e-л+h^ +

4+г

—д + ы+ ( — Л+

+

+ /4+г+1

д — Ы+^- e н+

9+ „—л+^+

Здесь для краткости записи использованы следующие обозначения:

2+ = 9 ((А+ — Л+) н+) , 9(*) = г/ (ехр(г) — 1) ,

д

Ы+

Т+

— N+

+

^ (ехр (—Л1 н+) — 1) + Л1

Ы+ =

1

Л+Л+'

N + =

Л+ + Л+

н+ (Л+Л+)

2 •

(12)

(13)

Выполняя аналогичные действия, получим диффузионный поток д на [хг—1, хг] в точке х = хг—0 (слева):

1 ^х

£=£¿-0

«г/Г ( Л- + |Г) — «г-1/Г( ^ eЛ-hг

+

где

4,г

д — ы— л- +

9.

н-

+ Т4,г-1

—д + Ы- ^ eгЛ2 н'

н-

Ы-

"н-

+ N-

Г (ехр (Л- н ) — ц — Л-

(14)

(15)

остальные обозначения введены выше.

В узле хг должны выполняться условия непрерывности функции м(х) и диффузионного потока д [2]:

д|ж=ж,-_0 д|ж=^+0

(16)

'£=£¿-0 '£=£¿+0 ' Ч1Ж=^_0

Приравнивая потоки (12) и (14) и проводя необходимые вычисления, получим каноническую форму разностных уравнений:

— СгМг + 6^+1 = — I = 1, N — 1,

/Г

н—9 e ;

Т4,г— 1

6 = н+ 9 e ;

д — ЫГ ^ eл-h

н-

Л Г , Т +

н— 9 + н+

9 + тг 9+ + т Л- — Л'Ч

+ т

4,г

—д + ы -

Л — + 9

Л2 + н—

+

с

Я

а

с

+г

-Я + М+(|+ - А1

+

+ /4+1+1

Я - м+^ е-л+ь+

Аппроксимацию краевых условий удобно получить на двухточечном шаблоне тем же способом, что и для внутренних узлов. При этом не будет нарушаться общая структура матрицы коэффициентов. Выделяя в (2) поток 5+ на [хо, х1] в точке х = хо и подставляя соответствующее выражение (12), получим сеточную аппроксимацию левого граничного условия (2) при г = 0:

-ОоПо + ЪоЩ = —¿о

(18)

51 (/

Ъо = 7+ 1- У е

/+ V ь

+

-л+ь 1

_ 51 Со = /+

+

/+ I Ь7 - А+

52,

¿о

4,о

-Я + М + ( - А+ ь1 1

+

+ /+

4,1

Я - М + е-л+ь1 ь1

5з-

Проводя аналогичные выкладки, получим двухточечную сеточную аппроксимацию правого граничного условия (3) при г = N:

амим-1 - смим

N,

(19)

ам

/1 V1 Ьм е

См

Р1 /1

/- А- +

ь

N

+ Р2,

¿м = / ^ /+ЛГ

-д + м- ( а- + ь-

N

+ -1

д - М- ^ ел2 Ьм

+ Рз-

Сложные выражения (17) - (19) несколько упрощаются, если считать источник кусочно-постоянным, /4 (х) « /- на (Хг-1, Хг) и /4 (х) ~ /+г на (хг, х^+1). Коэффициенты Яг, Ъг, Сг в этом случае остаются прежними, а принимает более простой вид

¿г

/

4,г

М

- (

ь-

1 - ел2

+ А-

+ /4+г

М+

+

Ь+

(1 - ел+ь+) - А+

Данная упрощенная схема применялась в [31] для расчета задач внешней аэродинамики.

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

Отметим некоторые свойства полученной схемы. Согласно (17), коэффициенты аг, Ъг, сг являются всегда положительными, так как у(г) > 0 для любых г. Кроме того, имеет место диагональное преобладание сг > аг + Ъг, г = , что обеспечивает устойчивость метода прогонки и монотонность разностной схемы [2]. В силу потокового способа построения она обладает свойством консервативности, поэтому обеспечивает строгое выполнение интегрального закона сохранения для исходного дифференциального уравнения (1). Принципиальным является следующее. После вычисления иг в узлах, используя формулу (11), при необходимости можно рассчитать распределение искомой функции и(х) и диффузионного потока 5(х) внутри сеточного интервала [хг, хг+1 ] и, следовательно, на всем отрезке [а,Ъ].

Таким образом, построенная разностная схема является точной для уравнения (1) в предположении произвольных кусочно-постоянных на [а, Ъ] коэффициентов /1(х) — /з(х) и кусочно-линейной правой части /4 (х). В случае более сложной их зависимости на сеточном интервале погрешность аппроксимации схемы составляет О(Ь2).

1

У

3. Асимптотика коэффициентов разностной схемы

Рассмотрим поведение коэффициентов разностной схемы (17) в случае малых значений определяющего сеточного параметра АЛ, А = шах(|А1|, |А2|). При А1Л ^ 0, А2Л ^ 0 имеем

#м =

и

и

и

еш -1

1 - 2 + 12 + о(и4) , и = (А2 - А1) л,

е ~

1 + и + ^ + О (и3) .

(20)

Удерживая необходимое по ходу анализа число членов в разложении (20), получим, что коэффициенты точной схемы (17) с погрешностью О (Л2) стремятся к следующим выражениям:

/1_ + / - Л_ _ /

Л- + /з 6 2

+

Ьг

/+ Л+

+/+ +Л+

6

12 2

, Л Л- ,+ Л+\

«г + Ьг - + /з^

^ = -

Л-

т

/4,г-1 + 2/4,г

3

3

+

л+ /

2

4 , г +1 . 2/4Тг

3

3

(21)

Наличие слагаемых /3Л/6 в формулах для коэффициентов аг, Ьг, а также вид выражения для указывают на принадлежность предельного вида точной схемы (17) к классу сплайновых аппроксимаций [33, 34]. В частности, при соответствующем выборе приближенных формул (7) для коэффициентов /1(х) — /з(ж) можно добиться тождественного совпадения со сплайновой схемой итерационно-интерполяционного метода [34]. Это означает, что пределом применимости сплайновых разностных схем [33, 34] является выполнение условия АЛ < 1, при котором справедливы разложения (20). Как следует из (21), полученная аппроксимация конвективного члена дифференциального уравнения соответствует центральным разностям. В этом нет ничего удивительного, так как в условиях существования разложений (20) АЛ < 1 выполняется известное ограничение на сеточное число Пекле Ре^ = ^Л// < 2.

Если в выражениях (20) ограничиться слагаемыми порядка и и считать источник кусочно-постоянным, то мы придем к традиционной разностной схеме конвективного тепломассообмена [35]:

/Г

Л-

2

Ьг

/+_ + /2 Л+ +2

«г + Ьг - /;

'з 2 + /з+ 2 )

^г = { 2 /4,г-1 + 2 /4+г+1

4. Анализ результатов численных расчетов

Для анализа точности предложенной схемы и ее сравнения с другими распространенными на практике сеточными аппроксимациями было рассмотрено несколько тестовых задач. Одна из них приведена в [30]:

^ж

1

И,е ^ж2

вт(пж), ж е (0,1) м(0) = м(1) = 0.

(22)

а

с

а

с

Таблица 1

Результаты решения задачи (22) на равномерной сетке

Номер узла г 1 2 3 4 5 6 7 8 9 10

Ие = 102 и ■ 104 157 559 1173 1950 2826 3731 4592 5338 5909 6259

А ■ 104 а 125 235 322 380 402 388 338 251 90 —590

б — 117 76 —209 260 —444 689 — 1053 1700 —2582 4132

в 99 189 262 312 335 330 296 237 158 64

г — 1 —4 —8 —13 —19 —25 —31 —36 —40 —42

Ие = 103 и ■ 104 132 511 1106 1870 2740 3646 4514 5275 5866 6240

А ■ 104 а 127 242 335 399 429 421 378 303 200 11

б —5030 233 —5253 499 —5489 794 —5747 1111 —6037 1446

г — 1 —3 —7 —13 —19 —25 —31 —36 —40 —42

Аналитическое решение уравнения (22) имеет вид Ие ч Ие2

и(х)

п2 + Ив2

8ш(пх) +

п (п2 + Ие2)

1 — ео8(пх) — 2

е-Яе(1-ж)_е-Яе

1е

-Яе

В табл. 1 и на рис. 3 приведены результаты численного решения задачи на равномерной сетке с шагом Ь = 1/11 (х0 = 0, х11 = 1) при Ие = 100 (Ие^ = 9.1) и Ие = 1000 (Ие^ = 91). В табл. 1 приняты следующие обозначения: и, и^ — точное и приближенное решения, А = (иЛ — и) — погрешность численного решения, строки "а", "б" относятся к аппроксимации конвективного слагаемого односторонними и центральными разностями, строка "в" — схеме Н. И. Булеева, Г. И. Тимохина [25] (совпадает со схемой А. М. Ильина [11]), сторока "г" — настоящей работе.

Рис. 3. Численное решение задачи (22) при различных значениях числа Ие: 1 — точное решение, 2 — центральные разности, 3 — противопотоковые разности, • — результаты расчета по схеме, предложенной авторами, Ь = 1/11; а — Ие^ = 9.1, б — Ие^ = 91.

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

предложенной схемы (строка "г") оказывается почти на порядок ниже одной из лучших аппроксимаций (строка "в") [11, 25]. Так как левая часть уравнения (22) описывается обеими схемами точно, улучшение погрешности численного решения напрямую связано с учетом линейной зависимости источника на сеточном интервале.

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

При фиксированном значении £ для последовательности сеток хг = гкк, г = 0,1,... , 1/кк, где кк = к/2к, к — шаг самого грубого разбиения, к = 0,1, 2,..., определяется максимум разности решений на двух сетках

Zk,e = max

к/2к к/2к+1 иг — и2г

(23)

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

Zk)£ < С (к/2к )Р , к = 0,1, 2,..., (24)

согласно общему принципу сходимости [13], обеспечивает равномерную сходимость с порядком р, если С и р не зависят от £. Из (24) следует, что при фиксированном £

Рк,е = 1с§2 ( Zk,e/ Zk+1)£) , (25)

а ^ ^к,£) является линейной функцией аргумента к.

Данная методика была реализована на следующей модельной задаче [13]:

£и'' + (1 + х2) и' — ((х — 0.5)2 + 2) и + 4 (3х2 — 3х + 1) ((х — 0.5)2 + 2) =0, (26)

и(0) = —1, и(1) = 0.

На рис. 4, а (сплошная кривая) хорошо видно, что решение уравнения (26) содержит пограничный слой толщиной порядка £ на левой границе. Несмотря на то, что шаг грубой сетки (к = 1/8) превышает его толщину (порядка 1/512), разностная схема практически точно воспроизводит в узлах решение и диффузионный поток (рис. 4 б). Использование локального решения (9) или (11) в качестве интерполирующей функции дает принципиальную возможность получить распределение искомой функции и диффузионного потока внутри сеточного интервала. На рис. 4 соответствующие значения в зоне пограничного слоя показаны крестиками. Как видно, имеет место хорошее совпадение с общим ходом точного решения. Расчеты на более мелких сетках еще сильнее уточняют результаты.

На рис. 5 показана зависимость ^ ^к£) от параметра измельчения сетки к. Для различных £ графики являются параллельными прямыми, что указывает на равномерную сходимость. В табл. 2 приведены результаты вычислений рк,£ в широком диапазоне изменения к и £. В последнем столбце дается средняя по всем к величина р£. Максимальное (при £ ~ 1) и минимальное (при £ ^ 0) значения р£ соответствуют порядку классической и равномерной сходимости. Анализ таблицы (последний столбец) показывает, что порядок р£ практически не зависит от величины £, что соответствует свойству равномерной сходимости разностной схемы по этому параметру.

Рис. 4. Графики функции п(х) и диффузионного потока еп'(х) в тестовой задаче (26): сплошная линия — точное решение, • — результаты расчета по схеме, предложенной авторами, х — интерполяция по формуле (9) в области погранслоя; Н = 1/8, е = 1/512.

1 ' 2 ' 3 ' 4 ' 5 ' 6 к

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

Таблица 2 Результаты решения модельной задачи (26)

к 0 1 2 3 4 Среднее

е Ре

1/2 2.004 2.001 2.001 2.000 2.000 2.00

1/4 1.996 1.996 2.001 1.999 2.000 2.00

1/8 1.978 2.000 1.989 2.050 1.957 1.99

1/16 1.993 2.001 2.000 2.002 1.989 2.00

1/32 1.939 2.000 1.999 1.998 2.010 1.99

1/64 1.912 1.992 1.999 1.996 2.006 1.98

1/128 1.924 1.966 2.022 1.992 2.006 1.98

1/256 1.913 1.975 1.996 2.029 2.001 1.98

1/512 1.905 1.961 2.004 2.013 2.038 1.98

Рис. 6. Графики функции и(х) (а) и диффузионного потока еи'(х) (б) в тестовой задаче (27): сплошная линия — точное решение, • — результаты расчета по схеме, предложенной авторами, Н = 0.1, е = 0.01.

В завершение рассмотрим модельную задачу [13], содержащую точку поворота, в которой происходит смена знака коэффициента при первой производной и в ней /2(х) = 0:

еи'' + 2хи' = 0, и(—1) = -1, и(1) = 2. (27)

Задача имеет точное решение

г

Ф(1Л/£) + 3Ф(хД^) . 2 Г _,2

и(х) = —5¥(Т77ё)— ■ =

е

Внутренний пограничный слой находится в точке поворота при х = 0 и имеет толщину порядка е1/2. Здесь решение резко изменяется от —1 до 2. Считаем, что точка поворота находится в одном из узлов сетки. Схема (17) - (19) дает правильные значения и(х) в узлах даже на грубой сетке при N = 2 —10. На рис. 6 приведены значения и(х) и диффузионного потока еи'(х) при е = 0.01 и N = 20 (к = 0.1). Возрастание числа узлов сетки N приводит к уменьшению различия численного и точного решений.

Заключение

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

Список литературы

[1] Теория тепломассообмена / Под ред. А. И. Леонтьева. Изд. 2-е. М.: Изд-во МГТУ им. Н. Э. Баумана, 1997. 683 с.

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

Марчук Г. И. Методы вычислительной математики. М.: Наука, 1980.

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

Вабищевич П. Н. Монотонные разностные схемы для задач конвекции-диффузии // Дифференц. уравнения. 1994. T. 30. C. 503-513.

Morton K. W. Numerical Solution of Convection-diffusion Problems. L.: Chapman Hall, 1996.

Ши Д. Численные методы в задачах теплообмена: Пер. с анг. М.: Мир, 1988. 544 с.

Патанкар С. Численные методы решения задач теплообмена и динамики жидкости М.: Энергоатомиздат, 1984. 152 с.

Флетчер К. Вычислительные методы в динамике жидкостей. Т. 1: Пер. с англ. М.: Мир, 1991. 504 с.

Васильева А. Б., Бутузов В. Ф. Асимптотическое разложение решений сингулярно возмущенных уравнений. М.: Наука, 1973.

Ильин А. М. Разностная схема для дифференциального уравнения с малым параметром при старшей производной // Мат. заметки. 1969. Т. 6, вып. 2. С. 237-248.

Бахвалов Н. С. К оптимизации методов решения краевых задач при наличии пограничного слоя // Журн. вычисл. математики и мат. физики. 1969. Т. 9, №4. С. 841-859.

Дулан Э., Миллер Дж., Шилдерс У. Равномерные численные методы решения задач с пограничным слоем. М.: Мир, 1983.

Шишкин Г. И. Сеточная аппроксимация сингулярно возмущенных эллиптических и параболических уравнений. Екатеринбург: Изд-во УО РАН, 1992.

Багаев Б.М., Шайдуров В. В. Сеточные методы решения задач с пограничным слоем. Новосибирск: Наука, СП РАН, 1998.

Roos H.G., Stynes M., Tobiska L. Numerical Methods for Singulary Pertubed Differntial Equation: Convection-diffusion and Flow Problem. Berlin, 1996.

Боглаев Ю. П. О численных методах решения сингулярно возмущенных задач // Дифференц. уравнения. 1985. Т. XXI, №10. С. 1804-1806.

Алексеевский М. В. Разностные схемы высокого порядка точности для сингулярно возмущенной краевой задачи // Дифференц. уравнения. 1981. Т. XVII, №7. С. 11711183.

Емельянов К. В. Усеченная разностная схема для линейной сингулярно возмущенной краевой задачи // Докл. АН СССР. 1982. Т. 282, №5. С. 1052-1055.

Сечин А. Ю. Численный метод высокого порядка точности для сингулярно возмущенной краевой задачи // Изв. вузов. Математика. 1983. №7. С. 75-80.

[21] Задорин А. И., Игнатьев В. Н. О численном решении уравнения с малым параметром при старшей производной // Журн. вычисл. математики и мат. физики. 1983. Т. 23, №3. С. 620-628.

[22] ЛисЕйкин В. Д. О численном решении сингулярно возмущенного уравнения с точкой поворота// Журн. вычисл. математики и мат. физики. 1984. Т. 24, №12. С. 18121826.

[23] Андреев В. Б., Савин И. А. О равномерной по малому параметру сходимости монотонной разностной схемы Самарского и ее модификации// Журн. вычисл. математики и мат. физики. 1995. Т. 35, №5. С. 739-752.

[24] Allen D. N., Southwell R. V. Relaxation methods applied to determine the motion, in to dimensions, of a viscous fluid past a fixed cylinder // Quart. J. Mech. and Appl. Math. 1955. Vol. VIII, pt 2. P. 129-145.

[25] Булеев Н. И., Тимухин Г. И. О численном решении уравнений гидродинамики для плоского потока вязкой несжимаемой жидкости // Изв. СО АН СССР. Сер. техн. наук. 1969. Вып. 1, №3. С. 14-24.

[26] Spalding D. B. A novel finite-difference formulation to differential expressions involving both first and second derivatives // Intern. J. Num. Methods Eng. 1972. Vol. 4. P. 55.

[27] Гущин В. А., ЩЕнников В. В. Об одной монотонной разностной схеме второго порядка точности // Журн. вычиол. математики и мат. физики. 1974. Т. 14, №3. С. 789-792.

[28] Эль-МистикАви Т. М., Верле М. Дж. Численный метод расчета пограничных слоев со вдувом — экспоненциальная разностная схема с ромбовидным шаблоном // Ракетная техника и космонавтика. 1978. Т. 16, №7. С. 138, 139.

[29] КАлис Х. Э. О применении некоторых монотонных разностных схем для решения эллиптических уравнений второго порядка // Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1985. Т. 16, №2. С. 65-80.

[30] Булеев Н. И. Пространственная модель турбулентного обмена. М.: Наука, 1989. 344 с.

[31] Белоусов В. Л., Головачев Ю.П., Земляков В. В. Неявная экспоненциальная разностная схема расчета сверхзвукового обтекания тел вязким газом // Мат. моделирование. 1994. Т. 6, №10. С. 66-76.

[32] Бронштейн И. Н., Семендяев К. А. Справочник по математике. М.: Наука, 1981. 718 с.

[33] Ильин В. П. О сплайновых решениях обыкновенных дифференциальных уравнений // Журн. вычиал. математики и мат. физики. Т. 18, №3. 1978. С. 620-627.

[34] Гришин А.М., Берцун В.Н., ЗинчЕнко В. И. Итерационно-интерполяционный метод и его приложения. Томск: Изд-во Томск. ун-та, 1981. 160 с.

[35] Дульнев Г. Н., Парфенов В. Г., СигАлов А. В. Применение ЭВМ для решения задач теплообмена. М.: Высшая школа, 1990.

Поступила в редакцию 5 октября 1999 г., в переработанном виде — 26 июля 2002 г.