CC BY
11
и
ШВ. И. Наац, Е. П. Ярцева
_____Численное исследование рекурсивных и итерационных алгоритмов..._
ННТЕМНГПКН
ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ РЕКУРСИВНЫХ И ИТЕРАЦИОННЫХ АЛГОРИТМОВ В ЗАДАЧЕ МОДЕЛИРОВАНИЯ ПЕРЕНОСА АЭРОЗОЛЕЙ В АТМОСФЕРЕ
В. И. Наац, Е. П. Ярцева
COMPUTATIONAL INVESTIGATION OF RECURSIVE AND ITERATIVE ALGORITHMS IN THE PROBLEM OF MODELING AEROSOL TRANSPORT IN THE ATMOSPHERE Naats V.I., Yartseva E.P.
The construction and computational investigation of recursive and iterative algorithms for modeling aerosol transport in the atmosphere are fulfilled; the comparative analysis of algorithms in computing experiment is made.
Key words: ecological monitoring, mathematical models of the transport theory, the concentration of contaminants.
Выполняется построение и численное исследование рекурсивных и итерационных алгоритмов для моделирования переноса аэрозолей в атмосфере, проводится сравнительный анализ алгоритмов в вычислительном эксперименте.
Ключевые слова: экологический мониторинг, математические модели теории переноса, концентрация загрязняющих примесей.
УДК 519.5, 519.6
Работа посвящена моделированию диффузного переноса загрязняющих аэрозольных примесей в пограничном слое атмосферы, построению соответствующих вычислительных алгоритмов и их численному исследованию. Подобный программно-вычислительный комплекс может быть использован в задачах оценки экологического состояния воздушной среды вблизи промышленных регионов. Перенос загрязняющих примесей (ЗВ) описывается так называемым параметризованным уравнением переноса - линейным дифференциальным уравнением в частных производных параболического типа [2] (одномерный вариант):
Э£( х, Г) -/-л а Э
- + а • q(x,t) + b • — (v(x,t) • q(x,t)) -
dt dX
. \
XS (x, t'), (1)
-e4
dx
K(x, t-) dq(x!) 4
Эх
^(х,0) = £ (х), д(о,Г )= £ (?), £(1,Г £ (Г).
(2)
в котором все поля (концентрация ЗВ £(Х, * ) > скорость ветра V (х, * ), турбулентность К (х, *) > представляемая коэффициентом турбулентной диффузии и источник примесей £(х, *)) нормированы, то есть значения х, Г е [0;1], £( х, ?) е[0;1],
V (х, Г )е[0;1], К (х, Г )е[0;1] и
£ (х,Г )е [0; 1]. В уравнении (1) нормировочные коэффициенты (параметры) имеют вид:
/ ч V • T K T S T
a(t) = a* T, b = max , в = max2 , £ = max . (3)
X X 2 ^max
Параметры (3) - безразмерные величины. Если заданы начальные и граничные условия для
q(x, t) (2) на W = [0;l]x [0;l] и заданы значения a , b, в, £ то решение уравнения (1) следу-
ет рассматривать как функцию этих параметров, а именно q(x, t,a, b,в,£) . Если V = Vmax и K = Kmax определены однозначно, то параметры Ь ив определяют X и T следующим образом:
X = Kmax Ь = K*-l T = Kmax b = K* • fi2
Vmax в V* в Vm2ax в V'2 в
Таким образом, определяются пространственно-временные интервалы в исходной модели (1) -(2). Кроме того, предполагается, что все исходные данные представлены массивами дискретных измерений, полученных в эксперименте с определенными погрешностями. Все это обусловливает необходимость построения эффективных вычислительных алгоритмов, обеспечивающих сходимость и устойчивость получаемых решений к погрешностям в исходных данных.
В настоящей работе выполняется построение и численное исследование двух алгоритмов. Первый построен на основе конечно - разностного метода (шеститочечная схема) [2, 3], а второй - на основе так называемого метода «интегральных уравнений» [1, 4]. Рассмотрим кратко данные алгоритмы.
1. Построение вычислительного алгоритма на основе конечно-разностного метода с использованием шеститочечной схемы включает в себя следующие шаги:
а) вводится сетка узлов {(х1:, tj)} на интервалах х е [0; 1], t е [0; 1], / = 0, т, ] = 0,п , х. = Ах- /, Ах = 1/т; t■ = At • ., At = 1/п.
I ■> J ;
б) вводятся обозначения д(х/, tj ) = д.., V(х/, tj) = ., К(х/, tj ) = К.., £(х/, tj ) = ,
1 = А^Ах2 и выполняется конечно-разностная аппроксимация производных на основе центральной разностной схемы и шеститочечной, входящих в уравнение [2]. В итоге вычислительный рекурсивный алгоритм получает следующее представление:
16 ( . ) 16 =
2- К/,.' +1+(1+16К.,.)- д.,.+1 2- К',] дг+и+1 =
= £Dt • Si, j +
+
к 1в г \ 1в г Л
—•Dx • Vi, j +—• (Ki _i,j - Ki+i, j)+—• Ki, j V У
'l_a • Dt _lb • Dx • (v +1 ,j _ V_i,j)_ 1в • Kj 4 ' lb-K,j+^K «,j _ K'_'j)_ lb 4
У
• qi_1, j +
•q,j +
+
• qi+i, j.
/ = 1, т-1, ] = 0, п-1. (4)
Вычислительная схема (4) представляет собой совокупность разностных уравнений, аппроксимирующих дифференциальное уравнение (1) во внутренних узлах сетки и является системой линейных алгебраических уравнений с трехдиагональной симметричной матрицей, для численного решения которой в работе применяется метод прогонки [2].
2. Вычислительный алгоритм на основе метода интегральных уравнений [1, 4] основан на предварительном интегрировании дифференциального уравнения (1) по переменной / и по-
следующему его сведению к интегро-дифференциальному уравнению. Подробный вывод и обоснование метода можно найти в работе [1]. Основные формулы вычислительного итерационного алгоритма представляют собой следующее:
Кк. = ехР |-£(0Г £•(а, •Ах2 + р•Ах• {ук,I- Vk-1,1))| =
= ехр<|- Я • £ Щ •Ах - (а1-Ах + р- (vk, г- ^-1, г ^ = Кк, г,.(1), (5)
Фк,. = дк ,0 • ехР |- Е ®г Ах2 •(а .■ •Ах 2+р •Ах •(^, .■- vk-1, .■ ^+
.=0
J At ~
+ Х Е А^Г £к, г'Кк, . •Ах 2 =
.=0 х
: дк ,0 • ехр]-1 • Е •(а. •Ах 2+р •Ах • (vk,.- vk-1,.))[+
. =0
.=0
'(V -1)
= -1Т (р • Ах • V*. - 6 • (Кк,. - Кк-,.)) • (дк:-1 - д^) )-
Ум Ах2
- 6 • К к. • (дк:;1! - 2дк:-1) + дк--1)) = У£Г’, (7)
Ах Ах
дй = Фк,.(1)-1 • Е Щ • Кк,.(1)• ¥к:-1). (8)
.=0
Условие выхода из итерационного вычислительного процесса формулируется так: вычисляется значение:
1 т п
р=—т—л ЕЕ1 С-. (9)
п • (т-1) к=^;‘=1| 1
и проверяется условие: если р <£, то д.. = д.(., в противном случае : = : +1.
3. Вычислительный эксперимент. Численное исследование описанных выше алгоритмов проводилось на основе тестового примера, подробно описанного в [3]. Условно данный тестовый пример называется «Блок исходных данных». В нем генерируются значения (массивы)
исходных данных {х. }т , {(,■ }п , {^. }тхп , {к.,. }тхп , {£.,. }тхп , а также точное решение д(х,t) ,
представленное массивом {д.. } ^ .
Точность получаемых приближенных решений каждого алгоритма определяется путем вычисления отклонения приближенного решения ~(х, t) от точного дТ (х, t):
1 т п .
а=—(—1) ХХд к- - ~к,
п • (т-1)
(10)
Вычисления проводились при следующих исходных данных V0 = 10 м/с,
К0 = 100 м/с2 , д0 = 0.001 кг/м3, £0 = 0.0001 кг/м"•с , а = 0.1, р = 0.1, 6 = 0.001, X = 1,
1 = 1/2, Ах = 0.05. В «Блоке исходных данных» вычислялись значения At = 0.00125, m = 20, n = 800, X = 972 м, T » 3.1 с, V » 31.5 м/с, K » 306 м/с2,
у у у у max / у max / у
q_ » 0.002 кг/м3, ^max » 0.000181 кг/м3-с.
Ниже приведены результаты расчетов по рекурсивному и итерационному алгоритмам, определяемых соотношениями (4) и (5) - (9), в виде массивов - массив q = {q. .}, i = 1, m, 2,
у = 1, и, 80, который представляет собой точное решение уравнения переноса:
г 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005
0.071 0.068 0.066 0.065 0.064 0.065 0.067 0.069 0.072 0.074 0.075
0.208 0.200 0.194 0.189 0.189 0.192 0.198 0.205 0.212 0.217 0.219
0.407 0.392 0.379 0.370 0.369 0.374 0.386 0.400 0.414 0.424 0.428
0.651 0.628 0.606 0.593 0.590 0.599 0.618 0.641 0.664 0.679 0.686
q = 0.924 0.890 0.860 0.841 0.837 0.850 0.876 0.909 0.941 0.964 0.972
1.203 1.159 1.121 1.096 1.091 1.108 1.142 1.185 1.226 1.256 1.267
1.469 1.416 1.368 1.338 1.332 1.352 1.394 1.446 1.497 1.533 1.547
1.702 1.639 1.585 1.549 1.543 1.566 1.614 1.675 1.734 1.776 1.791
1.883 1.815 1.753 1.714 1.707 1.733 1.786 1.854 1.918 1.965 1.982
1.999 1.927 1.862 1.820 1.812 1.840 1.897 1.968 2.037 2.087 2.105
массив q1 = kj}, j ,2 = = 1, n,80, полученный в результате расчетов, прово
рекурсивному алгоритму (первый алгоритм):
' 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005
0.071 0.068 0.066 0.064 0.064 0.065 0.067 0.069 0.072 0.074 0.075
0.208 0.199 0.193 0.189 0.188 0.191 0.197 0.204 0.212 0.217 0.219
0.407 0.389 0.377 0.369 0.367 0.373 0.385 0.399 0.413 0.424 0.428
0.651 0.620 0.605 0.592 0.589 0.598 0.616 0.639 0.662 0.681 0.686
~ = 0.927 0.878 0.859 0.840 0.836 0.849 0.874 0.907 0.939 0.968 0.972
1.203 1.129 1.119 1.096 1.091 1.106 1.139 1.181 1.223 1.264 1.267
1.469 1.368 1.366 1.339 1.333 1.351 1.391 1.443 1.494 1.548 1.547
1.702 1.575 1.579 1.552 1.544 1.565 1.612 1.672 1.731 1.798 1.791
1.883 1.735 1.744 1.716 1.708 1.733 1.785 1.852 1.918 1.995 1.982
1.999 1.837 1.848 1.822 1.812 1.841 1.897 1.968 2.039 2.123 2.105
и массив q2
= fe j} , i = 1, m,2, j = 1, n, 80
полученный в результате расчетов, проводимых по
итерационному алгоритму (второй алгоритм):
92 =
0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005
0.071 0.069 0.066 0.065 0.064 0.065 0.068 0.069 0.072 0.075 0.075
0.208 0.202 0.194 0.189 0.189 0.192 0.198 0.205 0.213 0.218 0.219
0.406 0.392 0.378 0.371 0.368 0.375 0.385 0.401 0.414 0.424 0.428
0.651 0.627 0.606 0.593 0.591 0.599 0.617 0.641 0.664 0.679 0.685
0.923 0.890 0.860 0.841 0.838 0.850 0.877 0.909 0.940 0.964 0.972
1.203 1.159 1.121 1.096 1.091 1.107 1.142 1.185 1.225 1.256 1.267
1.469 1.416 1.369 1.337 1.332 1.352 1.393 1.446 1.497 1.533 1.546
1.702 1.639 1.586 1.549 1.543 1.567 1.615 1.675 1.733 1.776 1.792
1.883 1.814 1.753 1.715 1.707 1.733 1.786 1.854 1.918 1.965 1.982
1.999 1.927 1.862 1.822 1.813 1.842 1.897 1.969 2.037 2.085 2.105
Погрешность решения о{дт, ~) (10) для первого алгоритма составила 0.02, для второго алгоритма - 0.003.
В вычислительном эксперименте исследовалась сходимость приближенных решений к точному решению для каждого вычислительного алгоритма. В ходе численных исследований было установлено, что важную роль в этом процессе играет соответствующий выбор значения параметра 1:
Таблица 1
Влияние параметра 1 на сходимость вычислительного процесса (о1 - погрешность для первого алгоритма, 0 2 - для второго алгоритма)
1 т п т х п о 1 0 2
1/2 10 200 2000 0.0296 0.0034
1/4 10 400 4000 0.0291 0.0035
1/6 10 600 6000 0.0289 0.0036
Расчеты показывают (таблица 1), что для повышения точности расчетов можно уменьшать 1, что в свою очередь приведет к резкому увеличению размерности задачи. Увеличивать размерность задачи нежелательно по двум причинам: а) поскольку исходные данные измерены с погрешностями, то большой объем данных может привести к неустойчивой работе алгоритмов; б) обеспечить модель большим объемом экспериментальных данных зачастую не представляется возможным. Поэтому желательно снизить объем вычислений, обеспечивая при этом приемлемую точность работы алгоритмов. В вычислительном эксперименте установлено, что минимальная погрешность решения, получаемого с помощью первого алгоритма, может быть достигнута порядка о » 0.01. При этом количество узлов составило 200 тыс. Дальнейшее увеличение числа узлов сетки приводит к увеличению погрешности за счет накопления погрешностей вычислений.
Для второго алгоритма (итерационного) минимальная погрешность решения составила 0.003 при числе узлов 2 тыс. В этом случае размерность задачи в 100 раз меньше, чем в предыдущем случае, а точность решения при этом на порядок выше.
Далее в вычислительном эксперименте исследовалась устойчивость получаемых решений к погрешностям в исходных данных. Алгоритм исследований на устойчивость подробно описан в [1]. Результаты соответствующих численных исследований представлены в таблице 2.
Таблица 2
Результаты исследования устойчивости вычислительных алгоритмов
Погрешность, вносимая в источник £ {х, ?) - правую часть исходного уравнения (1), е (%) Коэффициент усиления ошибки получаемого приближенного решения для первого алгоритма - ^ Коэффициент усиления ошибки получаемого приближенного решения для второго алгоритма - Т]2
0.1 0.03142 0.01228
2.5 0.78934 0.30684
5.0 1.58612 0.60296
7.5 2.38921 0.88877
10.0 3.19743 1.16555
Анализ результатов, представленных в таблице 2, свидетельствует о том, что оба алгоритма являются устойчивыми при определенных условиях. Так для первого алгоритма максимальное допустимое значение погрешности в исходных данных составляет до 2,5 %, а для метода интегральных уравнений - до 7,5 %, что говорит о большей устойчивости данного алгоритма. Значение коэффициента усиления ошибки больше единицы можно считать не приемлемым.
Таким образом, в результате численных исследований были выявлены особенности работы двух вычислительных алгоритмов, обосновано предпочтение алгоритма на основе метода интегральных уравнений по сравнению с алгоритмом на основе конечно - разностного метода. Эффективность итерационного алгоритма значительно выше, поскольку для его работы требуется значительно меньший объем вычислений и исходных данных, он обладает большей точностью получаемых решений и устойчивостью к погрешностям в исходных данных.
ЛИТЕРАТУРА
1. Наац В. И., Каргин Н. И. Итерационные методы численного решения задач переноса на основе интегральных уравнений // Известия вузов. Северо-Кавказский регион. Естественные науки. № 3. Приложение. Ростов-н/Д, 2004. - С. 3-16.
2. Наац В. И., Ярцева Е. П. Построение и исследование вычислительного алгоритма для параметризованного нестационарного уравнения переноса примесей в атмосфере // Современные достижения в науке и образовании: математика и информатика: материалы Международной научно-практи-
ческой конференции. — Архангельск, 2010. — С. 561—562.
3. Ярцева Е. П. Численные исследования алгоритма решения уравнения переноса субстанции в турбулентной атмосфере // Ломоносов — 2010: материалы XVII Международной конференции студентов, аспирантов и молодых ученых. — М.: МГУ, 2010.
4. Ярцева Е. П., Наац В. И. Итерационный метод численного решения уравнения переноса примесей в атмосфере // Актуальные задачи математического моделирования и информационных технологий: материалы VI Всероссийской открытой научно-практической конференции. — Сочи, 2010. — С. 178—180.
Об авторах
Наац Виктория Игоревна, ГОУ ВПО «Ставропольский государственный университет», доктор физико-математических наук, доцент, заведующая кафедрой математического анализа. Сфера научных интересов - математическое моделирование, численные методы, комплексы программ. УШаас @ yandex.ru
Ярцева Елена Павловна, ГОУ ВПО «Ставропольский государственный университет», доку-ментовед кафедры математического анализа. Сфера научных интересов - математическое моделирование, численные методы, комплексы программ.
уаЛ8еуа_е1епа @таЛ.ги
