CC BY
55
3. Merkulov V.I., Kharkov V.P., Shamarov N.N. Optimizatsiya kollektivnogo upravleniya gruppoy bespilotnykh letatel'nykh apparatov [Optimization of collective management of a group of unmanned aerial vehicles], Informatsionno-izmeritel'nye i upravlyayushchie sistemy [Information-measuring and control systems], 2012, No. 7, pp. 3-8.
4. Kharkov V.P., Merkulov V.I. Sintez algoritma ierarkhicheskogo upravleniya gruppoy bespilotnykh letatel'nykh apparatov [The synthesis algorithm of hierarchical control of a group of unmanned aerial vehicles], Informatsionno-izmeritel'nye i upravlyayushchie sistemy [Information-measuring and control systems], 2012, No. 8, pp. 61-67.
5. Gayduk A.R., Kapustyan S.G. Kontseptsiya postroeniya sistem kollektivnogo upravleniya bespilotnymi letatel'nymi apparatami [The concept of building systems for collective management of unmanned aerial vehicles], Informatsionno-izmeritel'nye i upravlyayushchie sistemy [Information-measuring and control systems], 2012, No. 7, pp. 8-16.
6. Roytenberg Ya.N. Avtomaticheskoe upravlenie [Automatic control]. Moscow: Nauka, 1992, 576 p.
7. Brayson A., Kho Yushi. Prikladnaya teoriya optimal'nogo upravleniya [Applied optimal control theory], Translation from English. Moscow: Mir, 1972, 544 p.
8. Chernous'ko F.A., Kolmanovskiy V.B. Optimal'noe upravlenie pri sluchaynykh vozmushcheniyakh [Optimal control with random perturbations]. Moscow: Nauka, 1978, 352 p.
9. Merkulov V.I., Drogalin V. V., Lepin V.N. i dr. Aviatsionnye sistemy radioupravleniya [Aircraft radio control system]. Vol. 1. Printsipy postroeniya sistem radioupravleniya. Osnovy sinteza i analiza [The principles of radio systems. Basics of synthesis and analysis]. Moscow: Radiotekhnika, 2003.
10. Merkulov V.I., Drogalin V.V. i dr. Aviatsionnye sistemy radioupravleniya [Aircraft radio control system]. Vol. 2. Radioelektronnye sistemy samonavedeniya [Electronic guidance system]. Moscow: Radiotekhnika, 2003.
Статью рекомендовал к опубликованию д.т.н., профессор В.П. Харьков.
Меркулов Владимир Иванович - ОАО «Концерн «Вега»; e-mail: [email protected]; 121170,
Москва, Кутузовский проспект, 34; тел.: 84992499476; д.т.н.; профессор; заместитель генерального конструктора.
Миляков Денис Александрович - к.т.н.; начальник лаборатории.
Самодов Игорь Олегович - тел.: 89653479164; инженер.
Merkulov Vladimir Ivanovich - JSC "Radio Engineering Corporation "Vega"; e-mail:
[email protected]; 34, Kutuzovskiy prospekt, Moscow, 121170, Russia; phone: +74992499476; dr.
of eng. sc.; professor; deputy of general constructor.
Milyakov Denis Alexandrovich - cand. of eng. sc.; head of laboratory.
Samodov Igor Olegovich - phone: +79653479164; engineer.
УДК 519.6:532.5
И.Н. Шабас, Л.Г. Чикина, А.Л. Чикин
ВЛИЯНИЕ ГРАНИЧНЫХ УСЛОВИЙ Ш-ГО РОДА НА УСТОЙЧИВОСТЬ ЗАДАЧ КОНВЕКЦИИ-ДИФФУЗИИ ПРИ ПРОТИВОПОТОКОВОЙ АППРОКСИМАЦИИ*
Представлены результаты разностной аппроксимации трехмерной задачи конвекции-диффузии с краевыми условиями третьего рода. Приводятся представления оператора конвекции-диффузии, приводящие к М-матричности. Формулируются условия устойчивости разностной схемы задачи конвекции-диффузии противопотоковой аппроксимации
*
Работа выполнена при финансовой поддержке Министерства образования и науки РФ (проект №1420, государственное задание вузов, базовая часть).
конвективных членов. Приводятся оценки решения на основе принципа максимума. За основу теоретического исследования взята методология математического моделирования и вычислительного эксперимента, предложенная академиком А.А. Самарским и развитая в работах российских и зарубежных исследователей. При аппроксимации задачи конечными разностями необходимо сохранить основные свойства исходных дифференциальных операторов. Поэтому при пространственной аппроксимации уравнения конвекции-диффузии, в котором конвективная часть записана в недивергентной форме, выбрана противопотоковая схема (ППС). Проведеное численное тестирование подтвердило выод о том, что наличие граничных условий третьего рода оказывает влияние на устойчивость в завасимости от величины х, входящей в это условие в виде коэффициента при свободном члене при фиксированном значении числа Пекле.
Задача конвекции-диффузии; противопотоковая разностная схема; граничные условия третьего рода; оценка решения.
I.N. Shabas, L.G. Chikina, A.L. Chikin
THE INFLUENCE OF THE THIRD KIND BOUNDARY CONDITIONS ON THE STABILITY OF CONVECTION-DIFFUSION PROBLEM WITH UPWIND APPROXIMATION
The results of the difference approximation of the three-dimensional convection-diffusion problem with boundary conditions of the third kind. Given representation of convection-diffusion operator, leading to M-matrix. We formulate conditions for the stability of the difference scheme convection-diffusion problem upwind approximation of the convective terms. Estimates of solutions based on the maximum principle. For a basis of theoretical research methodology derived mathematical modeling and computational experiment proposed by Academician A.A. Samara and developed in the works of Russian and foreign researchers. When approximating the problem by finite differences is necessary to preserve the main properties of the original differential operators. Therefore, when the spatial approximation of the convection-diffusion equation in which the con-vective part is recorded in non-divergence form, selected upwind scheme (US). Numerical testing. Numerical testing confirmed vyod that the presence of the boundary conditions of the third kind affects the stability depending on the coefficient of entering into these conditions in the form factor of the free term of a fixed value of the Peclet number.
Convection-diffusion problems; antiperspirant difference scheme; boundary conditions of the third kind; the assessment decision.
Введение. Широкое использование методов конечных разностей для решения задач математической физики вызвало необходимость детального изучения тех свойств разностных уравнений, которые непосредственно влияют на качество разностных схем. Такими свойствами, прежде всего, являются устойчивость и сходимость.
Схемы, которые в разностной форме сохраняют законы массы и энергии, называются консервативными [5]. Консервативные схемы дают более точные результаты. Свойство консервативности не связано напрямую с порядком точности схемы, кроме того, консервативные схемы первого порядка дают более точные результаты, чем неконсервативные схемы второго порядка [4]. Примером консервативности разностных схем являются противопотоковые разностные схемы [3, 4, 8]. В этих схемах используется односторонняя разность по пространству. При положительных скоростях используются левые разности (разности назад), а при отрицательных -правые разности (разности вперед по потоку). Противопотоковые разностные схемы вносят в решение схемную искусственную вязкость, и это следует учитывать при оценке точности результатов. К преимуществам схемы следует отнести свойство транспортивности [4]. Это свойство означает, что возмущения, возникающие в рассматриваемом процессе, передаются только в направлении скоростей за счет конвекции. Однако выбор разностей против потока отнюдь не
всегда гарантирует свойство транспортивности схемы [2, 4]. Что касается модифицированных противопотоковых схем [7], то такие схемы действительно точны, однако для слабо несамосопряженных задач.
Одна из серьезных трудностей, возникающих при решении уравнений конвекции-диффузии с преобладающими конвективными членами, заключается в возникновении пограничных слоев [4, 12]. Проблема погранслоя сыграла, вероятно, во многом решающую роль в развитии численных методов гидродинамики и, в частности, в решении вопроса о том, какими разностями - центральными или про-тивопотоковыми - следует аппроксимировать первые производные. С другой стороны, несмотря на первый порядок аппроксимации, разности «против потока» обладают многими полезными свойствами, такими как выполнение принципа максимума, монотонность, М-матричность. Все это определило большой интерес к разработке схем с разностями «против потока» [3, 4]. В отличие от центральных, противопотоковые разности не чувствительны к погранслою и, вообще говоря, к неадекватности получаемого решения точному.
Техника операторного подхода, предложенная в [6, 7] и развитая в работах [1, 9], позволяет успешно работать с несамосопряженными задачами. В этих работах получены неулучшаемые оценки устойчивости для широких классов двух- и трехслойных разностных схем. Тем не менее, вопросы общей теории устойчивости разностных схем для несамосопряженных задач нельзя считать разработанными полностью.
1. Постановка задачи. Уравнения, описывающие конвективно-диффузионный перенос, в несжимаемой среде (= 0) могут иметь различные эквивалентные
формы. В ограниченной трехмерной области О с границей Г будем рассматривать нестационарное уравнение конвекции-диффузии, когда конвективный перенос имеет недивергентную форму:
+ТV= г
+ Ь 1, (1) Ох[0, Т ], О = О^Г, О = {х = (х, у, г)}, где L - линейный оператор, действующий в гильбертовом пространстве H, 5 е У, / е Е, пространства Y, F с областями определения элементов О, О, Г
соответственно, где О = О^> Г, О = |х = (х,у,г)} . Оператор конвективно-диффузионного переноса Ь = Ьв + Ьс + Ьр состоит из операторов диффузионного переноса
df„ dS Л д
Ld = -—\ /лху —
дх у дх) ду у ду
( ясЛ Цху
д_ д1
конвективного переноса
дБ дБ дБ Ьс Б = и— + V — + w— (2)
2 дх ду дг
и свободного члена Ь^Б = /Б .
Задача (1) дополняется начальными данными £ (х,0) = V0 (х) и краевыми
условиями третьего рода на границе Г области О:
дБ
И— + = г, х е Г, I > 0. (3)
дп
Здесь 5 = 5(х) - концентрация рассматриваемого вещества; 1 - время;
Лх = Лх, ¡Лг = , цъ = ¡1* - коэффициенты турбулентной диффузии вещества;
V = (у1 , У2, У3) = (и, V, м>) - вектор скорости; коэффициент ( = /3(х, у, г, 1) > 0 описывает взаимодействие вещества со средой. Будем предполагать, что граница области определения решения гладкая, достаточную гладкость имеют функции Л = л(Х),х = х(Х),г = г(Х), заданные на границах, и решение задачи
5 = 5 ( Х) обладает достаточной гладкостью.
2. Дифференциально-разностная задача. Расчетная область произвольной формы помещается в прямоугольный параллелепипед О, затем вводится равномерная по всем направлениям разностная сетка Ок = Ок ^ Г^, векторным параметром Н = (Нх, Н^, Н), где к, Н, Н - соответствующие шаги сетки вдоль осей ОХ, ОУ, OZ. После проведения индексации ячеек определяется ОА - множество внутренних узлов сетки и Гh - множество граничных узлов. Аппроксимация задач (1)—(3) проводится в два этапа. Сначала эта задача аппроксимируется в области Ок х по пространственным переменным, а затем - по времени.
Аппроксимация краевых условий означает снос границы Г рассматриваемой нерегулярной области О на кусочно-линейную границу Г сеточной области Ок.
При аппроксимации граничных условий третьего рода правыми или левыми разностями используется идеология противопотоковых схем, когда выбор направления аппроксимации производной зависит от знака составляющей вектора скорости У(и , V, w), участвующей в граничных условиях. Разностные аналоги краевых условий третьего рода на семиточечном шаблоне приведены в табл. 1:
Таблица 1
Соответствие между осями и разностными аналогами краевых условий
На оси ЖЕ Лж + ХжНх£ ЛЖ £ = Г , Ле + ХЕНХ£ ЛЕ £ = г К "0]'к К"1]к Ж К "ыхк К'
На оси SN ¡5 +Х5Ну£ Л5 Л Л . г Н -ш Н £1к = '5, н '"ук Н у ~1к '' Ну Ну у у
На оси ВТ Лв +хА „ Лв „ .. , Лт + хА „ Лт „ _.. Н к"*1 =в к ^ т 2 2 2 2
В общем виде аппроксимация граничного условия третьего рода (3) будет
иметь вид е 5 + р £ = г , х, t еГхй, где £, £ - значения концентра© © © Н 1
ции вещества соответственно в некоторых граничном и приграничном узлах, = (л@ + X®Н& )/ Н0, з = ~Л© / Н& , Г© = Г©, где " ©" заменяется на
символы W, S, B, T, N E, если граница области приходится соответственно на О-1), (0-1), (к-1), 0+1), 0+1), (к,+1)-й узлы семиточечного шаблона, а = х,у, 2 соответствует осям ЖЕ, 5М, ВТ.
Аппроксимация свободного члена Ьр^ = /33и начального условия осуществляется точно: = и0 (х,0), х, г е ПА х {0}.
В результате задаче (1)-(3) поставлен в соответствие дифференциально-
д$
Разностный аналог + 1кяк = /к (х, г), х, г е^хЦ, где опеРатоР
дг
Ьк = Ьт + Ьс}г + Ь^ с разностными операторами диффузионного, конвективного переноса и разностным аналогом функции взаимодействия вещества, соответственно, и Ц = Ц ^{0} .
Исключив решение в граничных точках области и аппроксимировав производную по времени, можно перейти к неявной операторно-разностной схеме
-—-+кс=(х,г), $0=ио, п > о. (4)
т
Здесь оператор ЬИ - это оператор ЬИ, в который уже включены граничные
условия. Оператор ЬИ, соответственно, представим в виде = Ьт + ЬсИ + Ьрк .
В результате аппроксимации пространственных производных получена разностная схема на семиточечном шаблоне
П (0)ук$1-1ук + °(0)ук$у-1к + В(0)ук$ук-1 + ^(0)ук$ук + Т (0)ук$ ук+1 +
Np s + Fp s = fn
N (0)iiksii+1k + F(0)iiksi+1 ik J iik '
"П
(0)ук$ у+1к 1 Е(0)ук$1+1 ]к ^ ук '
коэффициенты которой для внутренних точек сеточной области имеют следующий вид:
X + У + 2 +
р чр р _ аук ™ук
^(0)ук=-н1 - Их ' ^=-%- ну' В(0)у=- н2 - К '
. + a
у
Г\Р — ' +1Jk 'Jk I 'J+1k iik I 'Jk+1 iik I I I I I I ■' I . I V"| о
D(0)ijk ~2 + Г2 + Г2 + +~r , , + Pijk'
( К К К hx hy hz j
z — y — x —
ТР _ aijk+1 wijk MP - aij+lk I vijk Т7Р __ai+ljk , Ujk
1(°)yk h2 + h ' N(0),jk h2 + h > F(0),jk h2 + h >
hz hz hy hy hx hx
„X
uxy + uxy uxy + uxy Uz + uz
_ Hi—1 jk + Hjk y _ Hj—1k + Hijk z _ *ijk—1 + rijk
где а.., =---—, ^ =—^-—, ^ =
* ук 2 ук 2 у 2
Приведенная конечно-разностная схема аппроксимирует исходное уравнение с погрешностью О(И) [3]. Включим решение в граничных точках области Цк на
р г
базе разностных краевых условий $ +--— $ = —— и получим разностную
& & & &
задачу для пространственного оператора: = , где Ьр - разностный
оператор в области , учитывающий вклады краевых условий, где
ТР? =Шр? -I- Яр V -I- Г)р V
!клк " укл1-\ук + £укау-\к + Вуклук-1 + Щуклук
+ Т1/к$дк+1 + +1к + мк , /п /ук 2 3® ® (0)ук ,
®=Ж ,Б,В,Т,Ы,Е ё®
® р = (1 - )® р)г]к., ® = Ж, £, В,Т, Ы, Е, 3 - символ Кронеккера.
Из построения видно, что вклады от граничных условий вошли в диагональ:
Щк = ВР(0)ук- X 3® ~Рт® Ру ' Щр)ук =- X ® (0)ук + Рук ■
) ®=Ж ,В,Т, N ,Е ё ® ®=Ж ,В,Т ,Ы ,Е
При конечно-разностной аппроксимации и естественном упорядочивании узлов рассматриваемой сеточной области Ок (а именно (1ДД),(1,1,2),(1Д,3),... ), состоящей из (Ых -1) х (Ы -1) х (Ыг -1) узлов, из уравнения конвекции -
диффузии получается СЛАУ, матрица которой имеет специальную семидиагональную структуру, причем диагональные элементы матрицы
линейного оператора Ьр положительны, а внедиагональные элементы
отрицательны.
3. Устойчивость разностной схемы задачи конвекции-диффузии. На сетке
ит = {гп = пт, о<гп <Т, т>0, п>0} = оти{0},
где = {гп = пт, т >0, п > 0}, рассмотрим двухслойную схему с весами
г,п+1 — сп _ _
-+ аЬр*п+ + (1 -а)Ьр8п = /п, (5)
т
где Ьр - оператор (разностный аналог уравнения конвекции - диффузии с
включенными граничными условиями), действующий в вещественном
конечномерном пространстве Н к со скалярным произведением (,), а а -
числовой параметр.
Для матрицы оператора Ьр в (5), которая имеет положительные
диагональные элементы и удовлетворяет условиям диагонального преобладания, получены [5] достаточные условия устойчивости на новой методической основе -с применением понятия логарифмической нормы оператора. Сначала рассматривается двухслойная разностная схема для системы однородных дифференциальных уравнений. Полученные общие условия конкретизируются на примере двумерного модельного уравнения конвекции-диффузии с граничными условиями первого рода. Условия
( „ \
Рк- X 3®
®=Ж ,В,Т ,Ы ,Е
® р > 0
®(0)ук > 0 (6)
1 +
&© ]
7 = 1,2,-1, } = \,2,...,и -\,к = \,2,...,И -1
' ' ' х ~ ' ' ' у ' ' ' ' г
являются достаточными [11] для того, чтобы матрица оператора Ь имела
положительные диагональные элементы и диагональное преобладание, а в силу неотрицательности внедиагональных элементов являлась М-матрицей.
Теорема 1. Пусть для оператора Ьр выполнено условие (6). Тогда
разностная схема (5) безусловно устойчива при а = 1 и условно устойчива при 0 < а <1 в , если
( ( ( „Л 1
т<
1
1 -а
max
1<i< N
V V
ßk - I
( p ^
1 I я p® 1 + s@—
0=W ,S ,B,T, N ,E
V
g 0
0 P
0 (0)ijk
0
> 0 •
JJ
При этом для разностного решения верна априорная оценка
И <1И Hill fk
V k=0
у s ü® 0 p
I s0 ®(0)ijk
0=W ,S ,B,T ,N,E g0
4. Оценки решения на основе принципа максимума. Используем методику, предложенную в [6]. На основе принципа максимума получим априорные оценки решения через правую часть уравнения, граничные условия и начальные данные. Для дальнейших исследований удобно записать это уравнение (4) в виде,
п+1
разрешенном относительно .у :
г\ \
.+ щр +Р |.п+1 = -ЖР .п+1 - £Р .п+1 - ВР .п+1 - ТР .п+1 -
+ Щук + Рук \.ук Ж (0)укЛг-1ук ° (0)укЛу-1к В(00)ук.ук-1 1 (0)укЛ ук+1
V
т
(7)
1
- Np sn+1 - Ep sn+1 + — sn + Jn
N (0)ijksij+1k E(0)ijk s i+1 jk T sjk T J jk ■
Обозначим через Р точку - центральную точку шаблона Ш(Р),
состоящего из совокупности семи точек рк, Ру+\к, Рук+г Ш'(Р) -
окрестность точки Р , т.е. все точки шаблона Ш'(Р) за исключением точки Р , в итоге Ш'(Р) - это шесть точек Р^к, Ру±\к, Р>*± 1- Тогда уравнение (7) можно записать в виде А(Р).(Р) = X В(Р, О).(О) + Е(Р), Р еП, где
Ое Ш (Р)
А(Р) > 0, В(Р, О), Е (Р) - сеточные функции, определенные для всех Р, О еП:
A( P) =1 + Dp + ß1]k,
1
B(Р,Q,) =,-Bp,-Ti ^ F(P)-fi.
Пусть D(P) = A(P) - X B(P, Q)-' в (7):
Q-ЕШ' (Р)
D(P) = A(P) - X B(P, Q) = PiJk.
Qe Ш'(Р))
Введем обозначения ¿[s] = A(P)s(P) - X B(P, Q)s(Q), P eQ. Наша
Qe Ш (Р)
цель - оценка в сеточной норме C PfP = max I f (p) I решения
PeQ
неоднородного уравнения
при условиях
L [ s] = F (P)
(8)
A(P) > 0, B(P, Q) > 0, D(P) > 0 для всех P, Q EÜ (9)
Линейный оператор L называется монотонным оператором, если из условия L[s(P)] > 0 для всех P EÜ следует, что s(P) > 0 для всех P e Q. Поэтому разностные схемы, удовлетворяющие при всех P EÜ условиям (9), называются монотонными разностными схемами, а для оператора L эти условия обеспечивают его монотонность [6], выполнение принципа максимума и корректность разностной задачи (8) в сеточной норме C .
Предположим, что сетка Q связная. Пусть Q - некоторое связное
* о
подмножество множества Q, а Q - дополнение Q до Ü. Введем нормы
P/P = max | f (P)|, PfP* = max | f (P)|. Для правой части задачи
C PeQ C *
PeQ
* о
f = f + f на Q решение задачи (8) представим в виде суммы
* о
s(P) = s(P) + s (P) :
s ( P) - решение задачи
s (P)
- решение задачи
L[s] = F, P EQ, L[s] = 0, P eQ,
о *
L[s] = 0, P eQ
[Щ = Р, Р еЦ.
Для нестационарных уравнений = Р(Р) оказывается целесообразным при оценке произвести "расслоение" сетки на множества меньшего числа
измерений и оценивать на них 5 (оценки на слое). Пусть окрестность Ш узла Р = Р(х, ?я+1), п > 0 состоит из узлов, принадлежащих двум временным слоям
г = гп и г = tn+1, так что Ш(Р(х, ?п+1))= Ш п+1+Ш "п, где Ш'п+1 - множество узлов Q = п +1) еШ(Р(х, ги+1 )), Ш п - множество узлов О = п) еШ(Р(х, 1)) . Тогда уравнение Ця] = Р(Р) можно записать в виде
А(Рп+х)з(Рп+1) = X В(Рп+1, О)8(О) + Ф(Рп+1),
Ое 0 п+1
где Ф(Рп+1)= X В(Рп+1, ®*(О) + Р (Рп+1). Рп+1 = Р(Х гп+1) . ОбознаЧим
Ое <
д'(Рп+1) = А(Рп+1) - X В(Рп+1,0). В дальнейшем понадобятся следующие
Ое 0 п+1
предположения:
д (Рп+1) > 0 VPn+1 е Ц А(Рп+1)> 0, В(Рп+1,0) > 0 УО еШ ^ УО е Ш "п, X в(Рп+1,0) > 0, X в(Рп+1,0) ^ 1 +
Ое Ш " д (Рп+1) Ое Ш"
пп
*
<
<
где c = const >0 не зависит от к,т. Предположим, что можно указать функцию
~ F(P )
F (P„+1) такую, что n+1
D\Pn+x)
< F(P , 1 где сеточная норма определена
II n+1 с,'
следующим обРазом: I\f(рп+1 Щ с = max |f(x, tn+i ^ = max |fn+i fk = f(x, tk).
h xe^ xeQ^
Выпишем коэффициенты уравнения (8) для узлов шаблона P :
о о * *
P e Q, n >0, P e Q, n = 0, P eQ, n > 0, P eQ, n = 0,
о *
P eQuQ, n = 0.
Для внутренней области Q при n >0 (8) примет вид
I — + Вр + В Ь^1 == -Шр чп+1 -БР чп+1 -БР чп+1 -тр чп+1 -
-МР сп+1 -ЕР 9п+1 ++1 + + /"п *Фш^+и + + т Ъ + +и.
п+1 °
Рассмотрим узел Р = из внутренней области Ц и его окрестность (а = 1,2,..7) при п >0
е=(сп оп+1 оп+1 оп+1 оп+1 оп+1 оп+1 \
а , -1Д , Яг+1 Д , 13 у-1к , Яу+1к , яук-1 , яук +1} '
Коэффициенты в (8):
А(Р) =—++Вук, д(Р) = А(Р)- X В(Р,О) = Вук,
Г Ое 0(Р)
В(Р О ) = ] 1 -жр -ЕР -яр -ыр -вр -Тр I
В(Р , ) 1 , "(0)ук, Е(0)ук, * (0)ук , * (0)ук, в(0)ук , Т (0) ук Г,
Р (Р) = УуШ.
Для внутренней области Ц при п = 0 (8) примет вид
Г1 + Вр +В ^к1 == -Шр V1 -*р V1 -Вр V1 -Тр V1 -
+ д(0)ук + Вук \яук ^(0)укяг-\ук * (0)укяу-1к в(0)укяук-1 Т(0)укЯук+1
)
-*р у1 - Ер V1 ++1 * + + 0
Окрестность узла Р = из внутренней области Ц составляет
О„(а = 1,2,..6): ={s)-lJk , , 4-1к , 4+1к , ^к-^ 4к+1}. КоэффиЦиенты в (8) в этом случае имеют вид
1
A(P) = -+D0) * ,
т
B(P Q ) = {-Fp -Sp -Np -Bp -Tp } B(P,ßi) { W(0)jk> F(0)ijk, S (0) ijk, N (0) jk, B(0)jk> T(0)jk
1 r, CV - S0 jk r0
D(P) = A(P)_ X B(P,ß) = - + , F(p) = -^ + j.
ße 0(P)
Для приграничной области П при п >0 (8) имеет вид
'1 ^ Р '
+ Щ(0)ук + Ж(0)ук + Вук
Vх ёж у
„^ __ _ пр „п+1 _ ТЗР ^п+1 _ ТР г.^1 _ 81]к °(00)ук8у-1к В(00)ук8ук-1 Т (00)ук 8¡у к+1
1 Г
-ыр 8п+1 -ЕР 8п+1 + +18п___Ж-Шр + /п
Ы (00)ук8у +1к Е(00)ук81+1ук ++ 81ук (00)ук + ¿ук-
т ёЖ
*
Узел Р = 8™+ из приграничной области П при п > 0 имеет окрестность
О (а = 1 2 7) : О ={8п 8п+1 8п+1 8п+1 8п+1 8п+1 8п+1 } °а (а 1,2,-- 7) : Оа {81ук,81-1 ук,81+1 ук , 8 у-1к , 8у+1к , 8 ук-1, 81ук+1 }-
Коэффициенты в (8):
1 Р г
А(Р)=-+Щ(0)рк +Рк, р (Р)=,
ёЖ ёЖ
В(Р О ) = 11 0 -ЕР -£Р -ЫР -ВР -ТР
В (Р, Оа ) 1 ,0, Е(00)ук, £ (00)ук , Ы (00)ук , В (00)ук , Т (00)ук
Р
Щ( Р) = -(1-)Ж(Р) к +Р к
ёЖ
*
Для приграничной области П при п = 0 (8) имеет вид
'1 Р„ '
^ + Щ(0)1ук + ' (0)ук + рук V ёЖ у
с1 == - <?Р с1 - Яр с1 - Тр с1 81ук °(00)ук8у- 1к В(00)ук8ук-1 Т (00)ук8ук+1
1
Г
-Ыр V1 -ЕР V1 ++1 £___—Жр + /0
Ы ^ +1к +1к + + т £ 0 у ё Ж(0) у + /к
Ж
При п = 0 окрестность узла Р = еП составляет О (а = 1,2,-.6):
0={?п+1 +1 \ а {81 -1 ук ,8 1+1 ук, 81у-1к, 8 ¡у+1к ,81ук-1, 81ук+1 } .
Коэффициенты в (8):
А(Р) =1+п(Р)к + Е^жрк +Вк, р(Р) = ^+у - ^-Щ^ук,
т ёж т ёж
Ш/(0)ук + рук,
Щ( Р) = 1-(1- ^ + В
т ё
° ж
В(Р, О1 ) = {0,-Е(0)ук ,-£(0)ук ,-Ы(0)ук ,-В(0)ук ,-Т(0)ук } .
Надо найти условия выполнения А(Р) > 0 и В(Р, О) > 0 с учетом, что В(Р, Оа) неотрицательны по построению для п > 0. Если показать справедливость неравенства Щ(Р) = А(Р)- X В(Р, О) > 0, то из
Ое Ш(Р)
неотрицательности коэффициентов В(Р, О ) будет следовать А(Р) > 0.
Рассмотрев случаи для внутренней области Ц и для приграничной области
*
Ц при п >0 и п = 0, для д(Р) и Р(Р) получим
D(P) =
ßk,
PeQ,
p
-1 TI
ßjk _(1_ ^ )W(0)jk, P eQ,
S„
PeQuQ,
n > 0, n > 0,
n = 0,
(10)
F (P) =
f
j'k '
f _■
ijk
s,
(0) jk >
P eQ,
*
P ei),
n > 0, n > 0,
(11)
S,
0ijk
P e QUQ, n = 0.
Представим в виде суммы правую часть Р (8) Р = Р + Р и решение этой
задачи s( P) = s (P) + s (P):
г (P) - решение задачи
P[s] = 0,P eQ
L[s] = F, P eQ,
(P) - решение задачи
* *
Z[s] = P, P gQ,
* о
¿[s] = 0, P gQ,
о о о
где P = У (P) для всех x GQ, n > 0 и p (P) = 0 в приграничных узлах *
X GQ, n > 0.
Используя методику [6], сведем задачу для s = s 0 + s 1 к задачам
* * * * * *
£[s0 ] = p0 (P), ] = p (P), где P0 и P1 удовлетворяют условию
* *
p (P) + p (P) = p(P). Запишем выражение для P0 и P1 : p (P) =
S
Р1 (Р) = I(Р) -^^ук, Ц, п > 0, Ц^Ц, п =
Теорема 2. Для решения нестационарного уравнения конвекции-диффузии (1), записанного в недивергентной форме, с краевыми условиями третьего рода в случае использования разностей "против потока" при выполнении условия
*
*
r
т
т
r
*
консервативности вещества ¡3 = 0 и условий (6) М-матричности оператора Ьр для неявной схемы (7) справедлива оценка
И, * (N1 + Т\\/II, +
тах
П>0
Л
1_ Е^
ё,
шр
ш
+
оот
ё,
шр
пО0)ук
ш у
. (12)
Доказательство. Условие М-матричности при 3 = 0 на О обеспечивает выполнение неравенства Б(Р) >0 по полученному условию (10):
* тах
Л
1_
Рш
ёш
шр
ш(0) т
шшр
ёш
о0)т
1 +
рш_ ёш
шр
ш(0) т
* Р^Р
(13)
где ||4 = тах |и(Р)|, ||£(
РеО
0г = шах\^0
к РеО,
^0 I,
к
С
с
*
+
И
И
0
к
С
с
С
С
Л
т
V ёш у
шр
ш(0) т
тах
*
Ре**
Л
(0) рк
V
V ёш у
шр
ш(0) рк
шшр „ ш (0)рк
&ш
(У шр ш(0) рк
V ёш
тах
РеП
шшР „ У (0)рк
&ш
[_ ршЛ (У шр ш(0) рк
V ёш
Неравенства (13) в сумме дают
* Р£п Р +
шах
ск п>0
Л
(
\
1_ РК
V ёш у
шр
ш (0)рк
+
шшр ш(0)рк
ёш
(
\
1_ рш
V ёш у
шр
П (0)рк
(14)
Чтобы оценить и , запишем уравнение Ь[И] = / (Р) иначе:
Ь[и] = И--Ь[и ] = Ф"+1, где Ф"+1 = — + /П. Так как Б =^>0, то
1
т
|ий+^ * ТФ
п+1
* И
+ т
НС:,
где
т
тах
|/(х)|. (15)
П
П
С
С
С
С
П
т
С
С
С
к
к
к
Суммирование (16) по n дает
, n+1
, откуда следует
ch k=0
* т
(16)
Суммируя (14) и (16), получим оценку (12).
Теорема 3. Для решения нестационарного уравнения конвекции-диффузии (1), записанного в недивергентной форме, с краевыми условиями третьего рода в случае использования разностей "против потока" при выполнении условия консервативности вещества /3 >0 и условий (6) М-матричности оператора Ьр для неявной схемы (7) справедлива оценка
INIc * PS0PC
+ тах
n>0
А
'ijk
r
+ max
о n>0
c
А
Vjk _
1_
\
pW
gW J
Wp
W (0)ijk
+
W
Sw
Wp
W(0)ijk
(
ßjk _
1_
Pw
S
Wp
W(0)jk
W J
(17)
Доказательство. Для области О О(Р) > 0, если /3 >0. Условие М*
матричности (6) при /3 >0 на О обеспечивает выполнение неравенства О(Р) > 0 по полученному (10). Таким образом, имеем выполнение О(Р) > 0 на всей области
О , следовательно из (10) и (11) получим (17).
5. Тестирование разностных схемы на модельных задачах. Влияние выбора различных разностных схем на сохранение свойства консервативности. В параллелепипеде построена сетка размером 21х 21x15. Шаг по горизонтали равен 2500 м, по вертикали - 1 м, шаг по времени 600 сек. Использованы четыре варианта поля скоростей: быстроменяющееся 1( sin2nx, —яуС082жх, —Ж2С0$2кх), однонаправленное 11(- 0.001,0,0), III - полученное из
гидродинамического расчета, IV - поле скоростей III, увеличенное в 10 раз. Для аппроксимации конвективных членов уравнения (2) применялась полностью неявная ЦРС, а для уравнения (1) - полностью неявная III 1С. Коэффициенты диффузии их и ¡и2 полагались либо равными константе, либо вычисленными по следующим формулам:
где R =
SP 0
du dz
+
0,00001, 0 < R < 0,25, 0,000001, R > 0,25, 0,001, R < 0,
k
N
c
h
N
C
h
*
c
C
— .
2
2
и Ху(х,у) = и0 Ь2 д/ 2(ды / дх)2 +(ды/ су + ду/ дх)2 + 2(ду / су )2 здесь g - ускорение силы тяжести; р0 - средняя плотность пресной воды; р - плотность соленой воды; и0 вычисляется по формуле и0 ~е11ъЬ41 3, где е полагался равным 5; и $ - числовая константа; Ь - характерный горизонтальный масштаб, равный шагу сетки.
Таблица 2
Влияние выбора различных разностных схем на сохранение свойства консервативности в (%)
Коэффициенты диффузии Поле скоростей ЦРС 1-я ППС 2-я ППС МППС
иху =1300 и2 =0,001 I 2,23 1,06 3,37 3,30
II 0,71 0,10 0,10 0,10
II 0,10 0,25 0,14 7,14
IV 15,511 6,17 16,98 52,32
иху =130 и2 =0,001 I 5,33 1,35 9,16 9,09
II 0,36 0,18 0,18 0,18
III 0,59 0,08 0,27 3,50
IV 45,23 44,48 40,46 60,76
Переменные по формулам I 5,70 4,33 7,18 6,92
II 0,65 0,09 0,09 0,09
III 1,37 1,25 1,45 6,69
IV 34,98 19,42 35,35 53,15
Проведенные тесты (табл. 1) показали, что центрально-разностная схема (ЦРС) дает меньшую потерю консервативности в случае вычисленного поля скоростей (III) при постоянных коэффициентах диффузии. Вторая противопотоковая схема (ППС-2) дает меньшую (по сравнению с остальными схемами) потерю консервативности при постоянных коэффициентах диффузии в случае вычисленного поля скоростей (IV), моделирующего большие значения вектора скорости. В остальных случаях среди рассмотренных схем меньшую потерю консервативности дает первая противопотоковая схема (ППС). Это позволяет сделать вывод о предпочтительности применения для расчетов первой противопотоковой схемы (ППС).
Влияние наличия граничного условия IIIрода на устойчивость. В единичном
квадрате О = [0,1]х[0,1] построена сетка размером 32 х 32. В качестве / иг из
уравнений (1) и (3) выбираются функции, удовлетворяющие аналитическому
решению £(', х) = в'$т(ш)$т(лу) задачи конвекции - диффузии
дБ 1 . „ дБ дБ
---ДБ + ы — + V— = /(х, у) .
д' Рв дх ду
При проведении расчетов были использованы четыре варианта поля скоростей: I(1. -1), П(1-2х, 2у-1), Ш(х+у, х-у), Щ$тлх, — ТГуС08Жх). Сходимость и ее скорость проверялась для двух случаев выбора граничных условий: с граничными условиями Ьго рода по всей области; со смешанными краевыми условиями: !-го рода по всем границам, кроме i — 1, на которой заданы граничные
условия Ш-го рода. В табл. 3-5 приведены: шаг по времени т, время счета t с, значения вычисления относительной погрешности (ОП) расчетов при применении ППС с граничными условиями I и III рода, изменение коэффициента х из граничного условия (3) на различных временных интервалах.
Таблица 3
Временной интервал [0; 20]
I рода I рода III рода (х=1) III рода (х=-29,65)
№ Pe т tc ОП % т tceK ОП % tc ОП % tceK ОП %
I 0,1 0,2 0,34 0,0283 20 0,03 0,4033 0,05 0,4567 0,03 111,1879
10 0,2 0,38 5,1307 20 0,04 15,8518 0,03 16,1592 0,03 15,839
1000 0,2 0,45 7,5169 20 0,04 24,2999 0,02 24,3009 0,03 24,2999
100000 0,2 0,49 7,5480 20 0,03 24,4221 0,04 24,4221 0,03 24,4221
II 0,1 0,2 0,28 0,0879 20 0,02 0,5195 0,03 0,5953 0,61 101,0861
10 0,2 0,34 0,7847 20 0,02 26,3385 0,02 27,1593 0,03 26,3086
1000 0,2 0,29 5,9855 20 0,02 63,0261 0,02 63,0266 0,04 63,0261
100000 0,2 0,28 9,2629 20 0,02 93,2668 0,03 93,2668 0,03 93,2668
III 0,1 0,2 0,4 0,0389 20 0,03 0,4674 0,03 0,5385 0,66 108,3246
10 0,2 0,38 2,4713 20 0,02 20,6347 0,03 21,2361 0,02 20,6095
1000 0,2 0,35 4,1846 20 0,02 37,2861 0,03 37,292 0,02 37,2859
100000 0,2 0,38 4,2163 20 0,03 37,7149 0,03 37,715 0,04 37,715
IV 0,1 0,2 0,41 0,0170 20 0,03 0,4457 0,03 0,5043 0,69 34,3994
10 0,2 0,49 3,6478 20 0,03 18,9338 0,04 19,1875 0,06 18,9203
1000 0,2 0,44 6,6848 20 0,04 33,9267 0,04 33,9278 0,04 33,9267
100000 0,2 0,43 6,7507 20 0,03 34,2803 0,04 34,2803 0,04 34,2803
Таблица 4
Временной интервал [0; 20], III рода, х=-29,6500, т = 2 0 , Р е = 1 О
задача I II III IV
tceK 0,03 0,02 0,03 0,02
ОП % 111,1879 107,3456 55,7502 19,5662
Нарушено диагональное преобладание
Таблица 5
Временной интервал [0; 2], т =0,02, III рода, задача I, Pe=10
X -5.0 -3,0 -2,5 -2,0 -1.0 -0.9 -0.8 -0,7 -0.6 -0.5 -0.4 -0.3 -0,3 -0,3 -0,3 -0,3 0.1 10 100
tceK 0,3 0,6 0.5 0.5 0.4 0,4 0,3 0,3 0,3 0,3 0,3 0,3 0,3 0.3 0.3 0.3 0.3 0.3 0.3
ОП% 7.1 7.1 >100 79,0 62,1 49.8 7.6 7.2 7.1
На величину относительной погрешности вычисления существенно влияет величина х, входящая в граничное условие третьего рода (3) в виде коэффициента при свободном члене при фиксированном значении числа Пекле (коэффициент
при производной в (3) равен обратной величине Ре). При Ре=0,1 нарушение диагонального преобладания (7) происходит при хе [-200,-1], но тесты показали, что счетная неустойчивость наступает при хе [-40,-10], особенно при х=-29,65 (табл. 3) При Ре=10 нарушение диагонального преобладания (6) происходит при хе [-100,-10], счетная неустойчивость наступает при хе[-2,5;-0,29] (табл. 5).
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Гулин А.В. Устойчивость разностных схем и операторные неравенства // Дифференциальные уравнения. - 1979. - Т. 15, № 12. - С. 2238-2250.
2. Ильин В.П. Методы конечных разностей и конечных объемов для эллиптических уравнений. - Новосибирск: Изд-во Ин-та математики, 2000. - 345 с.
3. Крукиер Л.А., Чикина Л.Г. Некоторые вопросы использования противопотоковых разностных схем при инженерных расчетах загрязнения в мелких водоемах // Инженерно-физический журнал. - 1998. - Т. 71, № 2. - С. 349-352.
4. Роуч П. Вычислительная гидродинамика. - М.: Мир, 1980. - 618 с.
5. Самарский А.А. Введение в теорию разностных схем. - М.: Наука, 1971. - 553 с.
6. Самарский А.А. Некоторые вопросы теории разностных схем // ЖВМиМФ. - 1966. - Т. 6, № 4. - С. 665-682.
7. Самарский А.А. Теория разностных схем. - М.: Наука, 1989. - 432 с.
8. Самарский А.А., Вабищевич П.Н. Численные методы решения задач конвекции-диффузии. - М.: Изд-во УРСС, 1999. - 248 с.
9. Самарский А.А., Гулин А.В. Устойчивость разностных схем. - М.: Наука, 1973. - 416 с.
10. Самарский А.А., Михайлов А.П. Математическое моделирование: Идеи. Методы. Примеры. - М.: Наука, Физматлит, 1997. - 320 с.
11. Чикина Л.Г., Шабас И.Н. Условия диссипативности и М-матричности разностного оператора конвекции-диффузии с граничными условиями третьего рода // Вычислительные технологии. - 2005. - Т. 10, № 6.
12. Gresho P.M., Lee R.L. Don't suppress the wiggles - they're telling you something! // Computers Fluids. - 1981. - Vol. 9. - P. 223-253.
REFERENCES
1. GulinA.V. Ustoychivost' raznostnykh skhem i operatornye neravenstva [Stability of difference schemes and operator inequalities], Differentsial'nye uravneniya [Differential equations].
1979, Vol. 15, No. 12, pp. 2238-2250.
2. Il'in V.P. Metody konechnykh raznostey i konechnykh ob"emov dlya ellipticheskikh uravneniy [The methods of finite differences and finite volumes for elliptic equations]. Novosibirsk: Izd-vo In-ta matematiki, 2000, 345 p.
3. Krukier L.A., Chikina L.G. Nekotorye voprosy ispol'zovaniya protivopotokovykh raznostnykh skhem pri inzhenernykh raschetakh zagryazneniya v melkikh vodoemakh [Some questions use protivopotokami difference schemes for engineering calculations contamination in shallow waters], Inzhenerno-fizicheskiy zhurnal [Engineering-physical journal], 1998, Vol. 71, No. 2, pp. 349-352.
4. Rouch P. Vychislitel'naya gidrodmамика [Computational fluid dynamics]. Мoscow: Мп,
1980, 618 p.
5. Samarskiy A.A. Vvedenie v teoriyu raznostnykh skhem [Introduction to the theory of difference schemes]. Moscow: Nauka, 1971, 553 p.
6. Samarskiy A.A. Nekotorye voprosy teorii raznostnykh skhem [Some problems of the theory of difference schemes], ZhVMiMF [Journal of Computational Mathematics and Mathematical Physicist], 1966, Vol. 6, No. 4, pp. 665-682.
7. Samarskiy A.A. Teoriya raznostnykh skhem [The theory of difference schemes]. Moscow: Nauka, 1989, 432 p.
8. Samarskiy A.A., Vabishchevich P.N. Chislennye metody resheniya zadach konvektsii-diffuzii [Numerical methods for solving convection-diffusion]. Moscow: Izd-vo URSS, 1999, 248 p.
9. Samarskiy A.A., Gulin A.V. Ustoychivost' raznostnykh skhem [Stability of difference schemes]. Moscow: Nauka, 1973, 416 p.
10. Samarskiy A.A., Mikhaylov A.P. Matematicheskoe modelirovanie: Idei. Metody. Primery [Mathematical modeling: Ideas. Methods. Examples]. Moscow: Nauka, Fizmatlit, 1997, 320 p.
11. Chikina L.G., Shabas I.N. Usloviya dissipativnosti i M-matrichnosti raznostnogo operatora konvektsii-diffuzii s granichnymi usloviyami tret'ego roda [Conditions dissipatively and M-metricnet differential operator convection-diffusion with boundary conditions of the third kind], Vychislitel'nye tekhnologii [Computational technologies], 2005, Vol. 10, No. 6.
12. Gresho P.M., Lee R.L. Don't suppress the wiggles - they're telling you something!, Computers Fluids, 1981, Vol. 9, pp. 223-253.
Статью рекомендовал к опубликованию д.ф.-м.н., профессор Ю.В. Тютюнов.
Шабас Ирина Николаевна - Южный федеральный университет; e-mail: [email protected]; 344090, г.Ростов-на-Дону, пр.Стачки, 200/1, корп.2, к.214; с.н.с.
Чикина Любовь Григорьевна - e-mail: [email protected]; кафедра высокопроизводительных вычислений и информационно-коммуникационных технологий факультета математики, механики и компьютерных наук; профессор.
Чикин Алексей Львович - Институт аридных зон Южного научного центра; e-mail: [email protected]; главный научный сотрудник.
Shabas Irina Nikolaevna - Southern Federal University; e-mail: [email protected]; 344090, Rostov-on-Don, st. Stachki, 200/1, b.2, a.214; senior researcher.
Chikina Lubov Grigoryevna - e-mail: [email protected]; the department VV IKT, faculty of Mathematics and Computer Science; professor.
Chikin Alexey L'vovich - Institute of Arid Ecosystems, South Science Center, Russian Academy; e-mail: [email protected]; chief scientific.
