Численное исследование сеточных моделей для нестационарного уравнения переноса Текст научной статьи по специальности «Математика»

УДК 517.2

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ СЕТОЧНЫХ МОДЕЛЕЙ ДЛЯ НЕСТАЦИОНАРНОГО УРАВНЕНИЯ ПЕРЕНОСА

© 2004 г. Н.И. Каргин, В.И. Наац

The construction of net models for the unsteadily-state equation of source transferring in the frontier layer of the atmosphere, which is based on explicit and implicit difference scheme, to carry out. The test example for realizing calculation experiment is developed. The results of numerical research of calculation algorithm convergence and stability are receiving.

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

Одномерное нестационарное уравнение переноса примесей в турбулентной атмосфере имеет вид [1]: дq(t, x) / дt + а(; )q(t, x) + д(у ((, x)q(t, x)) / дx --д(К((,x)дq(t,x)/дx)/дx = S((,x), (1)

q(t0, x) = р(x) , q(t, x0) = ф(() , q(t,X) = у/(() , t 6 [;о,Т] , X 6 [о,X], в котором q(t, х) - концентрация примесей; V(;, х), К (;, х) и S(;, х) - соответственно поля скорости ветра, турбулентной диффузии и источник примесей; а(;) - параметр, описывающий процессы вывода или привнесения примесей. Требуется построить вычислительную схему решения краевой задачи (1) относительно q(t, х) по известным исходным данным V(;, х), К(;, х), S(;, х), а(;), начальным и граничным условиям. В данной работе вычислительные модели строятся на основе явной и неявной разностных схем.

Свойства вычислительных алгоритмов сеточных моделей для первой краевой задачи уравнения теплопроводности с постоянными коэффициентами хорошо изучены [2]. Этого нельзя сказать о нестационарных уравнениях типа (1), в которые помимо исходных данных, представленных функциями V (;, х), К (;, х) и

S(;, х), а(;), входят также и первые производные V'x (;, х), К'х (;, х). В результате в вычислительный

алгоритм для уравнения (1) необходимо включить разностную аппроксимацию первой и второй производных искомого решения ¿¡; (;, х) , 4х (;, х), ^х (;, х) и разностную аппроксимацию первых производных исходных данных V'.,, (;, х), К'х (;, х). Все это усложняет

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

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

Параметризация уравнения переноса, построение для него сеточных моделей

Построение сеточных моделей для нестационарного одномерного уравнения переноса с начальными и граничными условиями (1) предваряется нормировкой всех переменных и функций, входящих в него. Способ нормирования уравнения (1) подробно описан в [1]. Далее (1) можно преобразовать к виду дq(t, х) / дх + а(;, х^(;, х) + Ь(;, х) дq(t, х)/ дх -

- с((, х) д2q(t, х)/дх2 = ё (;, х), (2)

где а(;,х) = {а(;) + в дV((,х)/дх};

Ь(;, х) = {в ■ V (;, х) -у дК (;, х)/дх};

с(;, х) = {у ■ К(;, х)}, ё(;,х) = £, ■ S(;,х);

нормировочные коэффициенты вычисляются по фор* 1 * о

мулам: а(() = Т ■«((); в = V ТХ ; у= К ТХ ;

* *

£ = S Т / <1 .

Все переменные и распределения, входящие в (2),

- нормированные величины, значения которых изменяются в диапазоне [0,1]. В итоге уравнение (2),

включая параметры а((), в, у и £, представляет собой параметризованный вариант задачи (1).

Построение явной разностной схемы решения уравнения (2) начинается с введения равномерной сетки узлов {(;у , х,)} на интервалах ;е [0,1] и

х 6 [0,1], у = 0,п , I = 0,т +1 и задания шаблона, т.е. множества узлов точек сетки, участвующих в аппроксимации дифференциального выражения. Шаг сетки:

Д; = п“1 и Д х = (т +1) 1. Внутренние узлы сетки: {(;у , х )}, ] = 1, п; I = 1, т, внешние узлы сетки:

{;0 = 0,х,}, I = 0,т +1; { , х0 = 0}, { ,хт+1 = 1},

] = 0, п . Чтобы аппроксимировать уравнение (2) в точке (;у , xi), введем шаблон, состоящий из четырех узлов {(;У, х, _! I , х,), , х,+! I {(у+ь х,)}, соответст-

вующий так называемой явной схеме сеточной модели [2]. Исходные распределения q(t, х), V(;, х), К(;, х), S(;, х), а(;) заменим приближенными значе-

1S

ниями сеточных функций q(tj, х,) = qji,

V ((у, х) = V, , К ((у, х) = К, , S ((у, х) = SJ1,

а ((у) = ау . Производные, входящие в (2), в точке ((у, х,) заменим разностными отношениями. Тогда

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

(д)-1(+и - qJ,,)+ а0Л + Ьу,(Дх)-1( - qу,,-1)-

i(к 2 J"1 (qj.i+l- 2qji + qu-1 )=dj

(З)

где у = 0, п -1, , = 1, т ; начальные распределения

q(tо = 0, х,) = qо , , , = 0, т +1; д(/у, х0 = 0) = Ц/ 0 и

q(tj■, хт+1 = 1) = qJ■,m+l, / = 0, п известны из начальных и граничных условий (1), а коэффициенты получают следующее представление:

а/, = а/ +в(Дх)-1 ( -V/-,/-1 ^

Ьц =?■ Vу,-г(Дх)-1(кл -Ки-1),

с/, = У К/, , =#■ SJ7 .

Выражение (3) преобразуем к виду:

qу+1,, = Бу, - (1 ■ Лу, - 1) ■ qу, +

+Я • Bji • q j,i-і + Я • Cji • q j,i+l

(4)

где Л/, = а у, ■Дх + Ьр ■Дх + 2 ■ с у,, ВJi = Ь у, ■ Дх - 2 ■ с у ,

С/, = су,, Б у, = ё/, ■ Д( , 1 = Д / Дх2 .

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

Неявная разностная схема для решения уравнения (2) определяется шаблоном, состоящим из четырех

узлов: (К , { ^ ((р. +1, { -1 ¡, +1, х, ^ (;J +1, х, +1 )} . В этом

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

q/+1,, - q j,i , qу+1,, - qу +1,,-1

■+а/+l,iqj+1,, + Ь/+,,-

дt

Дx

- cj+1,i (к2 )- (qj+1,i+1 - 2qj+1,i + qj+1,i-1 )- dj+1,i, (5) где aj +1,i =aj +1 +ß-(x)-1 (Vj +1,i - Vj +1,i-1 ) ,

bj+l,i = ß • Vj+l,i - r • (Дх)-1 (Kj+l,i - Kj+l,i-1 ), cj+1,i = у ' Kj+1,i , dj+1,i = % • Sj+1,i .

Выражение (5) преобразуем к виду Я Bj+1,i • qj+1,i-1- (1 + Я Aj+1,i )qj+1,i +

+Я • Cj+1,i • q j+1,i+1 = -(Dj+1,i + q ji ) !

2

Aj+li = aj+li ■ dx + b,-+i i -Д: + 2• c

j +l,i

cj+1, i

Bj+1,i = bj + 1,i ' Дx cj+1,i , Cj + 1,i = cj + 1,i :

(б)

(7)

Бу+1, = ё/+1, ■ Д(, 1 = Д / Дх , у = 0, п -1, , = 1, т , где qо,i, qj,о и qJ■,m+1 известны из начальных и граничных условий (1).

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

агУг-1 + ЬгУг + сгУг+1 = /г, (8)

где г = 1, т ; / = 0, п -1 и коэффициенты системы получают представление

аг = 1' В/+1,г, Ьг =-(1 + 1' л/+1,г ^ сг =1 ■С/+1,г, /г = -(Б/+1,г + qjг),

у0 = q/+1,0 , Ут+1 = q/ +1,т+1 .

Алгоритм метода прогонки включает в себя два этапа: прямую и обратную прогонки.

Прямая прогонка состоит в определении коэффициентов {Лг} и В}, I = 1, т

Al =-^, Bl = A Al =-

Bl =

£l

Vі bl

dl- alBl-l

cl

alAl-l + bl

l = 2, m ;

alAl-l + bl

d1 = f1 - a1y0 , d2 = f2 , • • d1 = f1, —.,

dm-1 = fm-1, dm = ./m - cmym+1,

обратная - в нахождении {уг}, г = т,1,-1 по формулам

Уш

dm - amB

шпш-1

amAm—1 + bm

ym-1 Am-1ym + Bm-1 , —,

Уг = Лгуг+1 + Вг, —., у1 = Л1.У2 + В1.

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

аг Ф 0 и Ьг ^ 0 , (9)

\Ь\ > |а^ + |с^ , г = 1,т . (10)

Выполнение условий (9) и (10) гарантирует существование и единственность решения системы (8) и возможность нахождения этого решения методом прогонки. Выполнение условий (9) определяется исходными данными, значения которых всегда положительны по их физическому смыслу. Условие (10) фактически означает справедливость неравенств

c

1+1 ■ |Л /+1,г | >1 ■ \В/+1,г | +1 ■ |С/+1,г |, (11)

в которых значения коэффициентов определяются соотношениями (7). Очевидно, что выполнение условий (11) будет обеспечено соответствующим выбором значения параметра 1 в диапазоне 0 < 1 < 1.

Постановка и реализация вычислительного эксперимента

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

Положим, что исходные данные уравнения (1) задаются функциями

q(t, х) = qо (1 + 0,9 ■ бш(5( + 3)) ■ (1 + 0,9 ■ бш(5х + 5)), V(;, х) = ^(1 + 0,75 ■ бш(2( +1,5)) ■ (1 + 0,8 ■ бш(7х +1,5)), К(;, х) = К0(1 + 0,75 ■ бш(5( + 5)) ■ (1 + 0,75 ■ бш(7х + 3)), а(() = а0 (1 + 0,75 ■ бш(5( + 5)), (12)

где qо , Уо , Ко, ао - заданные константы. Функциям (12) поставим в соответствие сеточные функции

qJi = q(tj, х,), V/, = V ((/, х,), К/, = к (;,, х,),

а/ =а($/). Таким образом, массивы исходных данных, кроме массива значений функции источника S= S((/, х,), определены. Выполним нормирование

значений исходных данных следующим образом:

= q/, / ^= V, / V *, = К, / К *, (13)

* * *

где q , V и К - максимальные элементы массивов к/,}, К,} и {к/, }, /'=0,п, , = 0, т +1. Значения а/ = а(( /) не нормируются, так как в соответствии с физическим смыслом значение параметра а0 должно задаваться в диапазоне 0 < ао < 1.

Для определения значений функции источника S/, = S((/, х,) необходимо решить уравнение (1) относительно S, подставив в него исходные данные. Требуется также вычислить значения производных

p(t, x) =

дq(t, x) дt

g(t, x) =

дq(t, x) дx

r(t, x) =

д q(t, x) дx 2

и ((, х) = дV((, х) / дх, w(t, х) = дК((, х) / дх , входящих в

исходное уравнение. Их находим, непосредственно дифференцируя (12). Далее р((, х), g(;, х), г((, х),

и((,х) и w(t,х) приближенно заменим сеточными

функциями р/, = р((/, х,), gJi = g((/, х,), г/, = г((/, х, ),

и/, = и((/, х,), Wji = w(tJ■, х,) и выполним их нормиро-

*

вание аналогично (13), учитывая при этом, что р , * * * * g , г , и и w - максимальные элементы массивов

{р/,}, {gJi}, {г/,}, {и/,} и {wJi}, р=о,п,

, = 0, т +1. Подставляя нормированные значения сеточных функций { р/,} { gJi } { г/,} { и/,} и { WJi }

у = 0, п , і = 0, да +1 в исходное уравнение (1), получим выражение €уі + ар€р + - суі£уі = , в

котором

¿Г ~ . и*ЧТ € и = Т т€ А *Т

*

р

У 1 -*Q

aji = “aj +~ *^ji , bji = ~ г Vji 2 * "'ji ■

X • p

= K r T = T S

cjl v2 * , ajl * Sjl .

X • p p *

X • p X • p

Вычислив значения dji, определяем далее

* /О * *

Б ,і = djip / Т, = Б уі / Б , где Б - максимальный

элемент массива }, / = 0, п , , = 0, т +1.

Таким образом, массивы значений всех исходных

I а ,

i aivtiiti uuuuju.u. iTiavvxiuui jiia iviinn uvva

данных: ( € i Í Xi {a, } {€j,l іVjl■}, j {#,,},

] = 0, п , і = 0, да +1 получены; вычислены значения

* * * * нормировочных констант д , V , К и Б . Теперь

можно приступить к реализации алгоритмов сеточных моделей, полагая неизвестными значения {д^-},

у = 1, п , і = 1, да при известных значениях { }> {€ ,

{а] }, {Ту-}, {$],}, {$}і}, у = 0,п , і = 0,да +1. Известными считаются начальные {¿€0-}, і = 0, да +1 и

граничные {ду0 } да+1}, У = 0,п условия.

Программная реализация вычислительных алгоритмов, соответствующих сеточным моделям, проводилась при следующих значениях исходных данных: X = 350 м, Т0 = 15 м/с, К0 = 40 м/с -2, а0 = 0,5 с -1, Дх = 0,1, Я = 0,5 , Д/ = 0,005 , п = 200 , да +1 = 10 . Далее приводятся результаты расчетов в виде массивов. Массив д = {¿у і}, у = 1,п,20 , і = 1,да представляет собой точное решение; массив ~(1) ={~(1)у,і},

у = 1, п,20 , і = 1, да получен в результате расчетов, проводимых по явной схеме сеточной модели; массив ~(2) = {~(2)у,г'}, у = 1,п,20 , і = 1,да получен в резуль-

0,112 0,230 0,367 0,489 0,567 0,581 0,529 0,422 0,287'

0,067 0,138 0,220 0,294 0,341 0,349 0,317 0,253 0,172

0,031 0,063 0,101 0,135 0,156 0,160 0,146 0,116 0,079

0,012 0,024 0,038 0,051 0,059 0,061 0,055 0,044 0,030

0,015 0,030 0,048 0,064 0,074 0,076 0,069 0,055 0,037

0,039 0,079 0,127 0,169 0,196 0,201 0,182 0,145 0,099

0,078 0,161 0,256 0,341 0,395 0,405 0,369 0,294 0,200

0,124 0,253 0,404 0,539 0,624 0,640 0,582 0,464 0,316

0,163 0,335 0,534 0,712 0,826 0,846 0,769 0,614 0,417

0,188 0,386 0,615 0,820 0,951 0,975 0,886 0,707 0,481 ^

~(1).

0,113 0,229 0,368 0,495 0,573 0,580 0,518 0,408 0,277 ^

0,068 0,137 0,221 0,296 0,342 0,345 0,307 0,241 0,164

0,031 0,063 0,101 0,135 0,155 0,155 0,137 0,108 0,074

0,012 0,024 0,038 0,051 0,058 0,058 0,051 0,040 0,027

0,015 0,030 0,047 0,063 0,071 0,070 0,062 0,049 0,034

0,039 0,078 0,126 0,166 0,187 0,185 0,162 0,128 0,088

0,078 0,158 0,254 0,335 0,378 0,372 0,326 0,257 0,178

0,123 0,250 0,401 0,530 0,598 0,591 0,520 0,410 0,284

0,162 0,332 0,531 0,704 0,799 0,794 0,703 0,555 0,383

0,187 0,383 0,613 0,815 0,930 0,933 0,831 0,658 0,451

0,113 0,229 0,368 0,495 0,574 0,583 0,523 0,412 0,279 ^

0,069 0,133 0,229 0,365 0,402 0,358 0,284 0,194 0,105

0,032 0,060 0,105 0,181 0,178 0,150 0,116 0,076 0,037

0,012 0,023 0,039 0,063 0,058 0,050 0,039 0,026 0,013

0,015 0,029 0,047 0,064 0,067 0,063 0,053 0,039 0,024

0,039 0,077 0,124 0,163 0,174 0,163 0,137 0,102 0,066

0,078 0,154 0,248 0,325 0,345 0,322 0,269 0,200 0,130

0,121 0,241 0,390 0,509 0,539 0,502 0,419 0,310 0,202

0,159 0,318 0,516 0,676 0,715 0,665 0,554 0,407 0,262

ч 0,181 0,363 0,595 0,788 0,835 0,775 0,642 0,466 0,292,

тате расчетов, проводимых по неявной схеме сеточной модели с применением метода прогонки.

На рис. 1, 2 показано графическое представление точного решения и решения, вычисленного по алгоритму явной разностной схемы. Точность расчетов при этом характеризуется величинами

ст(1) =ст(д,~ = q(1)), ст(2) =ст(д,~ = д(2)), опреде-

ляющими отклонение в среднем точного решения от

~(2) .

1

расчетных значений искомого а(д, ~) = — , где

а

m

І = 1

qji- qji

J=0

. Для рассматриваемого примера

точность вычислений составила ст(1) = 1,27Е - 02 и г(2) = 5,10Е - 02 . На рис. 3 показаны графики функ-

а

ций q, q1 = д(1) и q2 = д(2) для фиксированного момента времени.

Для сеточных моделей на основе явной и неявной разностных схем исследована сходимость методов с

Рис. 1. Графическое представление функции q = { qу ^}

Рис. 2. Графическое представление функции ~(2) = { ~j2)}

течением времени. На рис. 4 дано распределение ошибки СТу1 и СТу2) , ] = 1, п по временной оси. Результаты расчетов показывают, что с течением времени точность получаемых решений начинает снижаться, т.е. процесс медленно «расходится». При этом следует отметить, что сеточный метод на основе явной разностной схемы в этом смысле «работает» лучше метода на основе неявной разностной схемы. Очевидно, что значение Т не должно быть слишком большим, необходимо выполнение условия Т < X , например, X и Т • V).

В ходе проведения вычислительного эксперимента исследовалось влияние параметра Я на сходимость вычислительных алгоритмов. В табл. 1 приведены

q

Рис. 4. Отклонение искомых решений от точного sigma\ = СТр и sigma2 = , ] = 1, п с течением вре-

мени

значения параметра Я = At / Ах2 , размерность сетки узлов и значения величин ст(1) = ст(д,~ = q(1)) и

ст(2) = ст(д, ~ = q(2)). Результаты расчетов, представленные в табл. 1, показывают, что уменьшение значения параметра Я повышает точность расчетов, однако это приводит к увеличению размерности сетки узлов и в конечном итоге к накоплению вычислительных ошибок. Поэтому приемлемыми значениями Я можно считать Я = 1/4, Я = 1/6 при конкретных значениях исходных данных.

Изучалась также устойчивость методов к возможным погрешностям правой части уравнения (1). Устойчивость сеточных методов на основе явной и неявной разностных схем исследовалась следующим образом. Пусть вместо точного значения правой части

исходного уравнения Б(0) = {б(0)} известно его приближенное значение с погрешностью 8, т.е. имеем Б(8) = {б(88}. Величина погрешности 8 определяет-

Таблица 1

Влияние параметра Я на сходимость сеточных моделей

Я m +1 n Явная схема а(1) Неявная схема а(2)

1I2 10 200 1,27Е-02 5,1026Е-02

1I4 10 400 6,37Е-03 5,0941Е-02

1I6 10 600 4,24Е-03 5,09115Е-02

1I8 10 800 3,18Е-03 5,0902Е-02

1I10 10 1000 2,54Е-03 5,0846Е-02

ся по формуле f =

Рис. 3. Графики функций q , q1 = q(1) и q2 = q(2) для фиксированного момента времени € = 0,5

1

n - m

11 (S(f1 - S J>)2

Ji

J=li=1

n - m

J =li=l

где

S(f) = S (0)(l + ©£), © = { Ji} , J = 0, n , i = 0, m +1

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

отклонения получаемого решения ~(8) в ходе реализации вычислительного алгоритма от точного значе-

ния q

п =

(0)

nm

—II(S f >-S

n - m

J=li=1

nm

— ££(?!■? >- Æ ’)2

n- m

J=1=1

Пусть П1 , П2 - коэффициенты усиления ошибки,

характеризующие степень отклонения решения ~(8)

от точного значения q(0), получаемого по модели с явной и неявной разностными схемами.

В табл. 2 приводятся расчетные величины СТ1, СТ2, П1 и П2 при различных значениях 8 . Для всех рассмотренных случаев (строк табл. 2) СТ1 <8 и СТ2 <8 , что может свидетельствовать об устойчивости вычислительных методов. Сравнивая коэффициенты П1 и

Таблица 2

Влияние величины погрешности правой части уравнения (1) на коэффициент усиления ошибки при вычислении искомого решения по алгоритмам сеточных моделей

f Явная схема Неявная схема

а1 П1 а2 П2

0 1,2734Е-02 0 5,1015Е-02 0

2,7Е-02 1,2753Е-02 0,48 5,1060Е-02 0,12

5,9Е-02 1,2898Е-02 1,01 5,1168Е-02 0,26

8,6Е-02 1,2995Е-02 1,46 5,1305Е-02 0,38

1,2Е-01 1,3247Е-02 2,06 5,1363Е-02 0,55

1,5Е-01 1,3515Е-02 2,44 5,2079Е-02 0,66

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

Таким образом, на основе проведенных численных исследований можно заключить, что построенные сеточные модели для нестационарного уравнения переноса «работают» достаточно хорошо, обладают устойчивостью и сходимостью. Однако для достижения приемлемой точности требуется большое количество узлов сетки, определяемое соответствующим выбором параметра Я . Это вызвано не-

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

Литература

1. Семенчин Е.А., Наац В.И., Наац И.Э. Математическое моделирование нестационарного переноса примеси в пограничном слое атмосферы. М., 2003.

2. Самарский А.А., Гулин А.В. Численные методы: Учеб. пособие для вузов. М., 1989.

Северо-Кавказский государственный технический университет, г Ставрополь____________________________10 ноября 2003 г.