ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2024 Математика и механика № 90
Tomsk State University Journal of Mathematics and Mechanics
МАТЕМАТИКА
М ATHEMATICS
Научная статья
УДК 519.6 MSC: 65L04, 65L12
doi: 10.17223/19988621/90/1
Специальная разностная схема для решения жестких краевых задач конвективно-диффузионного переноса
Валентин Георгиевич Зверев
Томский государственный университет, Томск, Россия, [email protected]
Аннотация. Конвективно-диффузионное уравнение переноса лежит в основе описания широкого круга процессов в механике сплошных сред. Доминирование конвекции над диффузией, знакопеременность коэффициента при первой производной приводят к образованию локальных пограничных и внутренних переходных слоев с большими градиентами функции, что создает серьезные трудности при численном анализе задачи классическими разностными схемами. Традиционная аппроксимация первой производной центральными разностями при больших числах Пекле приводит к осцилляциям и нарушению монотонности численного решения. Чтобы избежать этого, требуется сильное уменьшение шага сетки в областях узких зон с большими градиентами. Использование односторонних разностей сильно размазывает искомое решение из-за схемной вязкости и приводит к потере точности. Практические потребности решения жестких краевых задач требуют разработки и применения вычислительных технологий, обеспечивающих монотонность, точность и экономичность численного анализа. В данной работе предложена новая специальная разностная схема для численного решения жесткого конвективно-диффузионного уравнения переноса. Доминирующий конвективный член исключен из явного рассмотрения путем перехода к самосопряженной форме уравнения, что позволяет использовать известные методы численной аппроксимации. Для построения разностного аналога дифференциального уравнения на трехточечном шаблоне используется метод контрольного объема. Полученная схема является монотонной, консервативной. На тестовых примерах показаны большие возможности предложенной разностной схемы при больших числах Пекле на грубых сетках при решении жестких краевых задач конвективно-диффузионного переноса. Ключевые слова: конвективно-диффузионный перенос, разностная схема, метод контрольного объема, трехточечный шаблон, монотонность решения
Благодарности: Исследование выполнено при поддержке проекта № 0721-2020-0032 государственного задания Минобрнауки России.
Для цитирования: Зверев В.Г. Специальная разностная схема для решения жестких краевых задач конвективно-диффузионного переноса // Вестник Томского государственного университета. Математика и механика. 2024. № 90. С. 5-17. doi: 10.17223/19988621/90/1
© В.Г. Зверев, 2024
Original article
A SPECIAL DIFFERENCE SCHEME FOR SOLVING STIFF BOUNDARY VALUE PROBLEMS OF CONVECTIVE-DIFFUSION TRANSFER
Valentin G. Zverev
Tomsk State University, Tomsk, Russian Federation, [email protected]
Abstract. The convective-diffusion transfer equation is often found in problems of hydromechanics and heat and mass transfer. The dominance of convection over diffusion and the change in sign of the coefficient at the first derivative lead to the formation of boundary and internal layers with high gradients of the function. This creates serious difficulties in numerical analysis of the problem using traditional difference schemes. The traditional method of approximating the first derivative using central differences at high Peclet numbers can lead to oscillations and violate the monotonicity of the numerical solution. To avoid this problem, it is necessary to significantly reduce the size of grid cells in narrow areas with large gradients of the unknown function. The use of one-sided differences significantly smears the desired solution, due to the viscosity of the scheme, and leads to loss of accuracy. The practical need to solve stiff boundary value problems requires the development and use of computational technologies that guarantee monotonicity, accuracy, and cost-effectiveness in numerical analysis. In this paper, a new special difference scheme is proposed for the numerical solution of a stiff equation of convective-diffusion transfer. The dominant convective term is eliminated from explicit consideration by transforming the equation into self-adjoined form, which permits the use of well-known numerical approximation techniques. The control volume method is used to construct a difference analogue of a differential equation on a three-point template. The resulting scheme is monotonic and conservative. The test examples show great possibilities of the proposed difference scheme for large Peclet numbers on coarse grids in solving stiff boundary value problems of convective diffusion transfer.
Keywords: convective-diffusion transfer, difference scheme, control volume method, three-point template, solution's monotonicity
Acknowledgments: This research was carried out with the support by state task (project No. 0721-2020-0032) of the Russian Ministry of Education and Science.
For citation: Zverev, V.G. (2024) A special difference scheme for solving stiff boundary value problems of convective-diffusion transfer. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika - Tomsk State University Journal of Mathematics and Mechanics. 90. pp. 5-17. doi: 10.17223/19988621/90/1
Введение
Математическое описание широкого круга процессов тепломассообмена опирается на конвективно-диффузионное уравнение переноса. Основным способом его решения являются численные методы. Доминирование конвекции над диффузией, знакопеременность коэффициента при первой производной приводят
к образованию локальных областей с большими градиентами функции - пограничных и внутренних переходных слоев, что создает серьезные трудности при решении краевой задачи классическими разностными схемами [1, 2]. Их причиной является различный масштаб процессов диффузии, конвекции, источников и стоков. Математически это приводит к возникновению в задаче малого параметра при старшей производной и порождает сингулярно возмущенную задачу.
Появление осцилляций и нарушение монотонности численного решения дифференциального уравнения переноса является одной из проблем, возникающей при аппроксимации первой производной уравнения центральными разностями.
Такая аппроксимация приводит к потере диагонального преобладания, несимметричности матрицы коэффициентов системы разностных уравнений и ухудшению ее свойств. Для обеспечения монотонности решения требуется сильное уменьшение шага сетки в области узких зон с большими градиентами решения - пограничных слоев, что ведет к большим затратам процессорного времени.
Используемый на практике переход к односторонним разностям сильно размазывает узкие зоны искомого решения из-за схемной вязкости и в конечном счете приводит к потере точности. При больших числах Пекле она вообще превышает физическую вязкость (коэффициент при второй производной) исходного уравнения.
Способу аппроксимации конвективного слагаемого в уравнении переноса в литературе уделяется большое внимание. В [3] для этой цели применяются локальные весовые интерполяционные кубические сплайны, что дает определенные преимущества перед широко используемыми монотонизированными аппроксимациями второго или третьего порядка. Метод сопряженных операторов предложен в [4] для построения разностных схем, являющихся точным дискретным аналогом исходной краевой задачи для обыкновенного дифференциального уравнения второго порядка.
В современных вычислительных технологиях для обеспечения монотонности решения применяют принцип регуляризации разностных схем [1]. Обзор исследований по проблеме численного решения сингулярно возмущенных краевых задач на основе применения сеточных методов и разностных схем высокого порядка имеется в [5-9] и не является целью данной работы.
Практические потребности решения жестких краевых задач конвективно-диффузионного переноса, возникающих перед исследователями, ставят задачу разработки и применения вычислительных технологий, обеспечивающих монотонность, точность и экономичность численного анализа [9].
В данной работе получена новая разностная схема для решения стационарного конвективно-диффузионного уравнения переноса. В основе ее построения лежат исключение на математическом уровне слагаемого с первой производной и применение метода контрольного объема [2] к дивергентной форме уравнения для получения разностного аналога дифференциального уравнения на трехточечном шаблоне.
Математическая постановка задачи
Рассмотрим одномерное стационарное конвективно-диффузионное уравнение переноса на отрезке [а, Ь] с краевыми условиями первого рода [6]
^(/¡(х) + /2(х)^ + /з(х)и = /,(X), (1)
ах ^ ах) ах
х = а: и = и0, (2)
х = Ь: и = ик.
Здесь u(x) - искомая функция; f1(x)—f4(x) - коэффициенты уравнения. Считаем, что выполняются обычные для корректной постановки задачи условия f1(x) > е > 0 и Тэ(х) < 0, причем значение е может быть сколь угодно малым. Из выражения (1) путем варьирования коэффициентов Т1(х)—/4(х) можно получить необходимую форму дифференциального уравнения. Выбор граничных условий первого рода не является принципиальным и принят лишь для упрощения изложения основных этапов построения разностной схемы на трехточечном шаблоне.
Методика построения разностной схемы
Введем на отрезке [а, Ь | произвольную неравномерную сетку □ = { а = х0 < ... < х1 <х1 < хм <...<хн = й}, /' = О...Ы, й,.+1 = х,.+1 -х,., /г,. = х,. -х,._!.
Из теории обыкновенных дифференциальных уравнений следует [10], что слагаемое с первой производной может быть исключено из явного рассмотрения. Используем для этой цели функцию ф(х), в результате получим преобразованный вид, эквивалентный исходному уравнению (1) [10, 11]:
(x ^
1 d (,_,., ч du
, ч , , Ф(x)fi(x)— I + /з(x)u = /4(x), ф(x) = exp ф(x) dx ^ dx
' /2 (x) fl( x)
dx
(3)
где х, - произвольная точка (узел) по координате х. Убедиться в правильности (3) можно путем непосредственной проверки.
Умножим все слагаемые (3) на ф(х) > 0, получим самосопряженную (дивергентную) форму уравнения (1):
а|ф(хШх)| + ф(х)/з(х)и = ф(х)/4(х) . (4)
ах^ ах)
Нетрудно видеть, что теперь в (4) при коэффициенте диффузии появился множитель с сильно меняющейся функцией ф(х), описывающей взаимодействие второй и первой производных. Формат уравнения (4) уже не содержит в явном виде исходной трудности и позволяет использовать известные методы численной аппроксимации.
Для получения разностного аналога дифференциального уравнения проинтегрируем (4) от Хг-1/2 до хг+1/2, где индексы г - 1/2 и г + 1/2 соответствуют серединам отрезков [х,ч, х,] и [х,, х1+1] соответственно, получим
аи
1—.
ах
Ji+1/2 - Ji-1/2 + |ф( x)f3(x)udx = j Ф( x)f4( x)dx= J =Ф( x)f1( x) • (5)
Балансовое выражение (5) является основой для получения различного вида разностных схем в зависимости от предположения о характере поведения коэффициентов Тэ(х), Т4(х) уравнения. Рассмотрим некоторых из них.
x
V i
x
x
-1/2
-1/2
Аппроксимируем функции ,/(х) и ф(х) в промежуточных узлах следующим образом:
т г и+1 — Ы1 т г и— и-1
+1/2 * Фг +1/^.71,1 +1/2 , , -1/2 * Фг-1/2Л,! -1/2 Г~
(6)
г+1
Ф^ = ехр П ^ ехр | ^ ^ 1 -+ ^
М. х) / х)
■ = ехр(г ), г =
•/1,1+1/22
Ф™ =ехрI! т "хГехрГл,-,.,2
/2,г-1/2 ^ I .. -ч - /2,1 -1/2^ , 1 = ехр(-г ), г = - ,
/1,г-1/2 2
/1 ,г
■/1,| + /1,1+1 /■ _ /2,г + /2,1+1
2 ' /2,''+1/2 = 2
Здесь г - сеточный параметр, имеющий смысл числа Пекле. Верхние индексы «+» и «-» относятся к интервалам справа и слева от узла , соответственно. Значения коэффициентов / ,-1/2 и /2, ¿1/2 определяются аналогично по формуле среднего арифметического от значений на концах сеточного отрезка. В случае сильно меняющегося коэффициента/1(х) согласно [2] целесообразно использовать формулу среднего гармонического.
Рассмотрим наиболее простой случай, когда /з(х) = /,, и /4(х) = /4,, - кусочно-постоянные на сеточном отрезке х,-1/2 < х < х,+1/2, взятые в ¿-м узле. Значение неизвестной и также возьмем в ¿-м узле для усиления его влияния в дискретном аналоге. В итоге остается интеграл от «быстрой» функции ф(х), который может быть взят точно. С учетом введенных обозначений (6) получим
"¡+1/2 -I -¡+1/2
Ш/ = ! ф(х)3х = ! ф(x)dx + ! ф(х)ёх г
1 - ехр(-г ) , ^! ехр(г+) -1
2 г - 2 г+
= ъ. Ш- + Ьм Ш+ 2 2
(7)
1 - ехр(-г )
ехр(г+) -1
у- (г -) =-^-^, (г +) =
г г
где у-, - сеточные функции аргумента г. Выражения для интегралов с источниками с учетом (7) имеют вид:
хг+1/2 хг хг+1/2
! ф(х)/ (х)^х = ! ф(х)/ (х)^х + ! ф(х)/ (х)^х *
хг-1/2 хг-1/2 хг (8)
2 2
¡+1/2 I ¡+1/2
|ф(х)/4(х^х = |ф(х)/4(хУх + |ф(х)/4(x)dx */4,1
Ь - Ь+1 +
Ш Ш+
2 2
= /4,,Ш, •
В итоге, подставляя (6), (8) в (5), получим коэффициенты канонического вида разностной схемы во внутренних узлах, которые зависят от сеточных функций
ф(г), у(г):
х
х
¡-1/2
х
х
х
-1/2
-1/2
aiui-i - ciui + biui+i = -di , i =1 ,N -1 , (9)
f1,i-1/2 , f1, i+1/2 , , Г 1 Г
a =Фг-1/2 -, bi =Фг+1/2 ~,-, С = a, + b, - f3 ,i Vi, d, = -f4,i Vi •
A,
A,,
Представляет интерес другой важный случай, когда Тэ(х) и Т4(х) - кусочно-линейные функции на сеточных отрезках [х,ч, х,] и [х,, х1+1], например
г . (/3,г+1 — /3,г) / N /• / N /• (/4,г+1 — /4,г) ✓ ч
/3 (х) = /3,1 +-^-~ (х - х X /4 (х) = /4,г +-^-~ (х - х ) .
A
i +1
A
i+1
Тогда интеграл от источникаТ4(х) принимает вид:
х!+1/2 х1 | ф(х)/4(х)ах = | ф(х)
г , (f4,i f4,i-1 ) , 4
f4,i +-;-(x - xi)
dx +
■ Ф j Ф(x)
f-1y Г(z " ) + кг
r . (f4,i+1 f4.i) 4
f4,i +-+-(x - x,)
Ai+1
A Л-(z -) + ^ Л+ (z + ) 2 2
dxг
exp(-z )(z +1)]
§-(z-) =11 -exp(-z )(z +1)|, §+(z+) =
2
(z - )2
+ f4,i+1-^ Г (z +), (10) exp( z +)(z +-1) +1]
(z + )2
Л" (2-) = (2~)(2~ ), Л+ (г+ ) =У+ (2+ )(2+ ) . Здесь п(2) - сеточные функции источника, играющие роль весовых множителей при периферийном и центральном узлах шаблона. Интеграл от линейного источника Тэ(х) аналогичен (10):
\ /з,, _ г (z - ) + /з,,
i+1/2
J ф(x)/ (x)udx г
A A
A (z- ) + A+i (z + )
+ /з,,+1 + (z + )j
(11)
Подставляя (10), (11) в (5), с учетом (2), (6) получим вид специальной разностной схемы во внутренних узлах:
а1и1-1 — СМ + Ьiui+1 = — а, г = 1, N — 1, (12)
/1,г —1/2 , /и+1/2
а = Фг —1 / 2-,-, Ьг = Фг +1 / Г
с, = a + bi - /-1-^- (z- ) + f3 i
A,+1
A (z - ) + Af Л+ (z + )
+ /з,+1^ Г (z + )
dt = 1/4, i-1y (z- ) + f4i
A (z -)+Af л+(z+)
+ /4,,+1^Г(z + •
Трехточечные выражения (12) образуют систему (Ы - 1) линейных алгебраических уравнений с трехдиагональной матрицей с неизвестными и1, ..., иЫ-1. Она решается прямым экономичным методом прогонки [12], требующим О(Ы) арифметических действий.
x
x
i-1/2
Отметим некоторые важные свойства специальной разностной схемы конвективно-диффузионного переноса. Согласно (12), коэффициенты щ, Ъ, являются всегда положительными, так как сеточная функция ф(г) > 0 для любых г. Кроме того, в матрице коэффициентов имеет место диагональное преобладание С ^ а + Ь, , = 1, N -1, что обеспечивает устойчивость метода прогонки и монотонность разностной схемы [9]. Также она обладает консервативностью и обеспечивает интегральный закон сохранения, справедливый для дифференциального уравнения (1).
Посредством коэффициентов щ, Ъ, учитывается решение однородной части уравнения (1). Основу их выражений составляет сеточная функция ф(г), описывающая взаимодействие конвекции и диффузии. Коэффициент ^ определяет неоднородную часть уравнения. Нетрудно видеть, что источник /4(х) при линейной зависимости от координаты х берется во всех узлах шаблона. Через сеточные функции п(г) осуществляется влияние решения однородной части уравнения на источник.
Асимптотика коэффициентов разностной схемы
Представляют интерес асимптотические выражения коэффициентов специальной схемы (9), (12) при малых значениях сеточного параметра г ^ 0, когда становится справедливым применение традиционных разностных схем. Разложение сеточных функций в ряд Тейлора в этом случае имеет вид:
Ф,- = ехр(-2") = 1- 0((2-)2), фт/2 = ехр(2 + ) = 1 + 2 ++ 0((2 + )2) , (13)
у- = 1 -ехр(-2-) = 1 -^ + 0((7-)2), = еХр(" +) -1 = 1 + 2+ + 0((2+ )2) ,
2
1 [1 - ехр(-2-)(2- +1)] _ 1 = 2 (2 ~ )2 = 2
1 [ехр(2+)(2+ -1) +1] 1
^ 2) = 2-(Т?-= 2
2
1_ ] 2 3
'1 2+ — +-
23
+ 0(( 2 -),
+ 0(( 2+ )2),
Ч~ (2- ) = У- (2- )(2 " ) =
(2+ ) = у+ (2+ )-Г (2+ ) =
3 _
4 3
[ 3 2+ — +-
4 3
+ 0((2- )2),
+ 0(( 2+ )2.
Учитывая (13), видим, что в предельном случае коэффициенты упрощенной схемы (9) стремятся к выражениям
(14)
а,ы,-1 - с,и, + Ь,и,+1 = ,
, = 1, N-1,
■/1,1-1/2 /2,,-1/2 . */!,/+1/2 , /2,,+1/2
а = —,---,—, ь = —,--,
Л,+1 2
с< = а, + ь, - /ъ,
2 2
4 =- /А,
К , К1 2 2
2
Нетрудно видеть, что (14) представляет собой традиционную схему с центральными разностями для решения конвективно-диффузионного уравнения, записанную на неравномерной сетке.
Для схемы (12) в предельном случае имеем
аи—1 — Сиг + Ьи+1 = —аг, г = 1N — 1 (15)
/1,1—1/2 /2,г—1/2 ; /1,1+1/2 /2,1+1/2
а = —:---:—, Ь = —:--+ -
A,+1 2
, , Л /з,,-1 , з/з,,
с = а. + b ——,--I---
2 4 4
2 2
d =-J ^Ai^ + з/U
2 4 4
a, , A,+1 2 2
! Ai+1 /з,г +1 2 4
Ai+1 f4,i+1
~2 Г"
Как видно, здесь также присутствует аппроксимация первой производной центральными разностями, однако выражения для источников имеют более сложный вид, характерный для схем сплайновой аппроксимации [13].
Результаты расчетов и их анализ
Для анализа точности предложенной схемы (12) и сравнения с другими разностными аппроксимациями было рассмотрено несколько тестовых задач. Одна из них приведена в [11]:
du 1 d 2u
dx Re dx'
u(0) = u(1) = 0
и имеет аналитическое решение
= sin nx, x e (0,1), (16)
Re Re2 L . . Ae~Re(1-x) -e"Re)
u(x)= —--sin(nx) +---—U-cos(nx)-2 , ,
n2 + Re2 n(n2 + Re2) [ (1 -e~Re) J
которое показано сплошной кривой 1 на рис. 1. Увеличение числа Re при старшей производной приводит к возникновению области резкого изменения функции - пограничного слоя (ПС) на правом конце расчетной области.
В табл. 1 и на рис. 1 приведены результаты численного решения задачи на равномерной грубой сетке с шагом h = 1/11 (х0 = 0, хц = 1) при числах Re = 102 и Re = 103 , что соответствует сеточным значениям Reh = 9.1 и Reh = 91 соответственно. В обоих случаях ПС является подсеточным масштабом, и в его область не попадает ни одного расчетного узла.
В табл. 1 приняты следующие обозначения: u, uh - точное и приближенное решения, Ah = uh - u - погрешность численного решения. Строки а, б относятся к аппроксимации конвективного слагаемого односторонними и центральными разностями, в - к схеме Н.И. Булеева, Г.И. Тимухина [14] (совпадает со схемой А.М. Ильина [15]), г - к схеме (12) данной работы. Результаты численного решения указаны значками на рис. 1.
Из рис. 1 видно, что односторонние разности а плохо обрабатывают область погранслойного изменения функции на правом конце, схема b с центральными разностями, как и следует теоретически, приводит к сильным осцилляциям
+
решения. Ситуация с традиционной схемой становится еще более драматичной с увеличением числа Re.
Рис. 1. Численное решение задачи (16) при различных значениях числа Re: 102 (a), 103 (b): 1 - точное решение, 2 - центральные, 3 - противопотоковые разности, о - специальная схема (12), h = 1/11 Fig. 1. Numerical solution of problem (16) for different values of the Re number: Re = (a) 102 and (b) 103. (1) exact solution, (2) central, (3) upwind differences, о is the special scheme (12), h = 1/11
Таблица 1
Результаты численного решения задачи (16) различными схемами на грубой сетке
Номе р узла i 1 2 3 4 5 6 7 8 9 10
Re = 102 u • 104 157 559 1 173 1 950 2 826 3 731 4 592 5 338 5 909 6 259
А104 а 125 235 322 380 402 388 338 251 90 -590
б -117 76 -209 260 -444 689 1 053 1 700 -2 582 4 132
в 99 189 262 312 335 330 296 237 158 64
г -1 -3 -6 -11 -16 -21 -27 -31 -35 -37
Re = 103 u • 104 132 511 1 106 1 870 2 740 3 646 4 514 5 275 5 866 6 240
А104 а 127 242 335 399 429 421 378 303 200 11
б 5 030 233 -5 253 499 5 489 794 5 747 1 111 -6 037 1 446
г -1 -3 -7 -13 -19 -25 -31 -36 -40 -42
Ошибка предложенной схемы (строка г табл. 1) оказывается почти на порядок ниже одной их лучших специальных аппроксимаций в [11, 12] и по уровню погрешности соответствует результатам [6, 13, 14]. Улучшение точности численного решения напрямую связано с учетом линейной зависимости источника на сеточном интервале. Предложенная схема правильно воспроизводит решение при любом значении Яе.
Практический интерес представляет случай переменных коэффициентов дифференциального уравнения конвективно-диффузионного переноса. Для его анализа была рассмотрена следующая задача [5, 9, 16, 17]:
еи" + (1 + х2)и' - ((х - 0.5)2 + 2)ы + 4(3х2 - 3х +1)((х - 0.5)2 + 2) = 0 , (17) ы(0) =-1, ы(1) = 0 .
На рис. 2 (сплошная кривая) хорошо видно, что решение уравнения (17) содержит пограничный слой толщиной ~ е на левой границе. Несмотря на то, что шаг грубой сетки превышает его толщину ~ е = 1/512, разностная схема (12) практически точно воспроизводит в узлах решение как при И = 1/16, так и при И = 1/4.
2.01.51.00.50.0-0.5-1.0-
u(x)
О
2.01.51.00.50.0-0.5-1.0-
u(x)
О
0.4 0.6 a
0.4 0.6 b
1.0 x
Рис. 2. Расчет функции u(x) в задаче (17) при различном шаге сетки: h = 1/16 (a), h = 1/4 (b): сплошная линия - точное решение, о - специальная схема (12), 8 = 1/512 Fig. 2. Calculation of the function u(x) in problem (17) at a different grid step: h = (a) 1/16 and (b) 1/4. The solid line is the exact solution, о is the special scheme (12), 8 = 0.01
В общем случае коэффициент f2(x) при конвективном слагаемом в уравнении (1) может менять знак. Знакопеременность f2(x) делает возможным появление внутри расчетной области локальных областей с большими градиентами искомой функции в виде внутренних переходных слоев. Поэтому представляет практический интерес анализ возможностей предложенной специальной разностной схемы (12) конвективно-диффузионного переноса и в этом случае. Рассмотрим модельную задачу [5], содержащую точку поворота, в которой происходит смена знака коэффициента при первой производной, и в ней f2(x) = 0 [9, 16, 17]:
8u" + 2xu' = 0, x £ (-1,1) , u(-1) = -1, u(1) = 2 .
(18)
2 z
Ф(z) = f
лЫ '
dt.
(19)
Задача имеет точное решение
, ч Ф(1/л/е) + 3Ф(х/^е)
и( х) =--¡=-,
2Ф(1/л/Ю V» 0
Внутренний пограничный слой находится в окрестности точки поворота при х = 0 и имеет толщину ~ . Здесь решение резко изменяется от -1 до 2. Считаем, что точка поворота находится в одном из узлов сетки. Результаты расчетов показали, что схема (12) дает правильные значения и(х) в узлах даже на грубой равномерной сетке при малом количестве узлов N. На рис. 3 в качестве примера приведены значения и(х) при числе узлов N = 20 (И = 0.1) и N = 4 (И = 0.5). Малый параметр при старшей производной равняется е = 0.01. Видно, что хорошо вычисляется функция внутри каждого из пограничных слоев, а также во внешней области. Увеличение числа узлов сетки N приводит к уменьшению различия численного и точного решений.
0.0
0.2
0.8
.0
0.0
0.2
0.8
x
2.01.51.0 0.5 0.0 -0.5-1.0-
u(x)
0.0 a
2.01.51.00.50.0-0.5-1.0-
u(x)
0.0 b
Рис. 3. Расчет функции u(x) в задаче (18) с внутренней точкой поворота при различном шаге сетки: h = 0.1 (a), h = 0.5 (b): сплошная линия - точное решение, о - специальная схема (12),
8 = 0.01
Fig. 3. Calculation of the function u(x) in problem (18) with an internal pivot point at a different grid step: h = (a) 0.1 and (b) 0.5. A solid line is an exact solution, o is the special scheme (12),
8 = 0.01
-1.0
0.5
0.5
.0 x
.0
0.5
0.5
1.0 x
Таким образом, приведенные результаты тестовых расчетов подтверждают изложенный теоретический анализ и показывают большие возможности предложенной специальной разностной схемы при решении жестких краевых задач по сравнению с традиционными аналогами, использующими центральные и односторонние разности при аппроксимации первой производной уравнения.
Заключение
Предложена новая специальная разностная схема решения жестких краевых задач конвективно-диффузионного переноса для граничных условий первого рода. Исследована асимптотика ее коэффициентов при малом значении сеточного параметра и получена связь с известными в литературе разностными аппроксимациями. Схема является монотонной, консервативной, обеспечивает устойчивое получение численного решения при больших числах Пекле на грубых сетках и имеет хорошие перспективы дальнейшего развития.
Литература
1. Самарский А.А., Вабищевич П.Н. Численные методы решения задач конвекции-диффузии.
М. : Эдиториал УРСС, 1999.
2. Патанкар С. Численные методы решения задач тепломассообмена и динамики жидкости.
М. : Энергоатомиздат, 1984.
3. Семёнова А.А., Старченко А.В. Разностная схема для нестационарного уравнения перено-
са, построенная с использованием локальных весовых интерполяционных кубических сплайнов // Вестник Томского государственного университета. Математика и механика. 2017. № 49. С. 61-74. doi: 10.17223/19988621/49/6
4. Воеводин А. Ф. Метод сопряженных операторов для решения краевых задач для обыкно-
венных дифференциальных уравнений второго порядка // Сибирский журнал вычислительной математики. 2012. Т. 15, № 3. С. 250-260.
5. Дулан Э., Миллер Дж., Шилдерс У. Равномерные численные методы решения задач с по-
граничным слоем. М. : Мир, 1983.
6. Багаев Б.М., Карепова Е.Д., Шайдуров В.В. Сеточные методы решения задач с погра-
ничным слоем. Новосибирск : Наука, 2001. Ч. 2. 224 с.
7. Задорин А.И. Разностные схемы для задач с пограничным слоем. Омск : ОмГУ, 2002.
118 с.
8. Liseikin V.D., Karasulji'c S., Paasonen V.I. Numerical Grids and High-Order Schemes for
Problems with Boundary and Interior Layers. Novosibirsk : IPC Novosibirsk State University, 2021.
9. Зверев В.Г. Численные методы решения задач с пограничным слоем. Новосибирск :
Наука, 2017.
10. Камке Э. Справочник по обыкновенным дифференциальным уравнениям. М. : Наука, 1976.
11. БулеевН.И. Пространственная модель турбулентного обмена. М. : Наука, 1989.
12. Самарский А.А., НиколаевЕ.С. Методы решения сеточных уравнений. М. : Наука, 1978.
13. Гришин А.М., Берцун В.Н., Зинченко В.И. Итерационно-интерполяционный метод и его приложения. Томск : Изд-во Том. ун-та, 1981.
14. Булеев Н.И., Тимухин Г.И. О численном решении уравнений гидродинамики для плоского потока вязкой несжимаемой жидкости // Известия СО АН СССР. Сер. техн. наук. 1969. Вып. 1, № 3. С. 14-24.
15. Ильин А.М. Разностная схема для дифференциального уравнения с малым параметром при старшей производной // Математичекие заметки. 1969. Т. 6, вып. 2. С. 237-248.
16. Зверев В.Г., Гольдин В.Д. Разностная схема для решения конвективно-диффузионных задач тепломассообмена // Вычислительные технологии. 2002. Т. 7, № 6. С. 24-37.
17. Зверев В.Г. Об одной специальной разностной схеме для решения краевых задач тепломассообмена // Журнал вычислительной математики и математической физики. 2003. Т. 43, № 2. С. 265-278.
References
1. Samarskiy A.A., Vabishevich P.N. (1999) Chislennyye metody resheniya zadach konvektsii-
diffuzii [Numerical methods for solving convection-diffusion problems]. Moscow: URSS.
2. Patankar S. (1980) Numerical Heat Transfer and Fluid Flow. New York: Hemisphere.
3. Semenova A.A., Starchenko A.V. (2017) Raznostnaya skhema dlya nestatsionarnogo
uravneniya perenosa, postroyennaya s ispol'zovaniyem lokal'nykh vesovykh interpolyatsionnykh kubicheskikh splaynov [The finite-difference scheme for the unsteady convection-diffusion equation based on weighted local cubic spline interpolation]. Vestnik Tomskogo gosudarstven-nogo universiteta. Matematika i mekhanika - Tomsk State University Journal of Mathematics and Mechanics. 49. pp. 61-74. DOI: 10.17223/19988621/49/6.
4. Voevodin A.F. (2012) A method of adjoint operators for solving boundary value problems for
second-order ordinary different equations. Numerical Analysis and Applications. 5(3). pp. 204-212. DOI: 10.1134/S1995423912030020.
5. Doolan E., Miller J., Shilders W. (1980) Uniform Numerical Methods for Problems with
Initial and Boundary Layers. Dublin: Boole Press.
6. Bagaev B.M., Karepova E.D., Shaidurov V.V. (2001) Setochnyye metody resheniya zadach
s pogranichnym sloyem [Grid methods for solving problems with a boundary layer]. Part 2. Novosibirsk: Nauka.
7. Zadorin A.I. (2002) Raznostnyye skhemy dlya zadach s pogranichnym sloyem [Difference
schemes for boundary layer problems]. Omsk: Omsk State University.
8. Liseikin V.D., Karasuljic S., Paasonen V.I. (2021) Numerical Grids and High-Order Schemes for Problems with Boundary and Interior Layers. Novosibirsk: IPC Novosibirsk State University.
9. Zverev V.G. (2017) Chislennyye metody resheniya zadach spogranichnym sloyem [Numerical
methods for solving problems with a boundary layer]. Novosibirsk: Nauka.
10. Kamke E. (1976) Spravochnik po obyknovennym differentsial'nym uravneniyam [Handbook of ordinary differential equations]. Moscow: Nauka.
11. Buleev N.I. (1989) Prostranstvennaya model' turbulentnogo obmena [Spatial model of turbulent exchange]. Moscow: Nauka.
12. Samarskiy A.A., Nikolaev E.S. (1978) Metody resheniya setochnykh uravneniy [Methods for solving grid equations]. Moscow: Nauka.
13. Grishin A.M., Bertsun V.N., Zinchenko V.I. (1981) Iteratsionno-interpolyatsionnyy metod iyegoprilozheniya [Iterative-interpolation method and its applications]. Tomsk: Tomsk State University.
14. Buleev N.I., Timukhin G.I. (1969) O chislennom reshenii uravneniy gidrodinamiki dlya ploskogo potoka vyazkoy neszhimayemoy zhidkosti [On the numerical solution of the equations of hydrodynamics for a plane flow of a viscous incompressible fluid]. Izvestiya SO AN SSSR, Series of Technical Sciences. 1(3). pp. 14-24.
15. Il'in A.M. (1969) Differencing scheme for a differential equation with a small parameter affecting the highest derivative. Mathematical Notes. 6(2). pp. 596-602.
16. Zverev V.G., Gol'din V.D. (2002) Raznostnaya skhema dlya resheniya konvektivno-diffuzionnykh zadach teplomassoobmena [Finite-difference scheme for solving convection-diffusion problems of heat-mass exchange]. Vychislitel'nyye tekhnologii. 7(6). pp. 24-37.
17. Zverev V.G. (2003) A special difference scheme for boundary value problems in heat and mass transfer // Computational Mathematics and Mathematical Physics. 43(2), pp. 255-267.
Сведения об авторе:
Зверев Валентин Георгиевич - кандидат физико-математических наук, заведующий лабораторией Научно-исследовательского института прикладной математики и механики
Томского государственного университета (Томск, Россия). E-mail: [email protected]
Information about the author:
Zverev Valentin G. (Candidate of Physics and Mathematics, Head of laboratory, Research
Institute of Applied Mathematics and Mechanics of Tomsk State University, Tomsk, Russian
Federation). E-mail: [email protected]
The article was submitted 22.12.2023; accepted for publication 05.08.2024 Статья поступила в редакцию 22.12.2023; принята к публикации 05.08.2024