Методика расчета вертикальной составляющей коэффициента турбулентной диффузии Текст научной статьи по специальности «Математика»

УДК 519.644

МЕТОДИКА РАСЧЕТА ВЕРТИКАЛЬНОЙ СОСТАВЛЯЮЩЕЙ КОЭФФИЦИЕНТА ТУРБУЛЕНТНОЙ ДИФФУЗИИ

Семенчин Евгений Андреевич д. ф.-м. н., профессор, зав. кафедрой

Кузякина Марина Викторовна аспирант

Кафедра высшей алгебры и геометрии, Кубанский государственный университет, Краснодар, Россия

Предложена методика расчета вертикальной составляющей коэффициента турбулентной диффузии в математической модели рассеяния примеси в приземном слое атмосферы

Ключевые слова: ФИЛЬТРАЦИЯ, КОНЦЕНТРАЦИЯ ПРИМЕСИ, ВЕРТИКАЛЬНАЯ СОСТАВЛЯЮЩАЯ КОЭФФИЦИЕНТА ТУРБУЛЕНТНОЙ ДИФФУЗИИ

UDC 519.644

A TECHNIQUE FOR COMPUTING OF THE TURBULENT DIFFUSION COEFFICIENT VERTICAL COMPONENT

Semenchin Evgeny Andreyevich

Dr. Sci. (Phys.-Math.), professor, Head of department

Kuzyakina Marina Viktorovna postgraduate student

The higher algebra and geometry department, Kuban State University, Krasnodar, Russia

The technique for computing of the turbulent diffusion coefficient vertical component in the context of a mathematical model of admixture dispersion in the surface layer is proposed

Keywords: FILTRATION, ADMIXTURE CONCENTRATION, THE TURBULENT DIFFUSION COEFFICIENT VERTICAL COMPONENT

Введение

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

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

Математическая модель, описывающая процесс рассеяния примеси в приземном слое турбулентной атмосферы имеет вид [1]:

где q(t, х, у, г) — средняя концентрация примеси в точке (х, у, г) е Е+ ,

Е+ ={( х, у, г) : х, у е (-¥; ¥), 2 > 0}, в момент времени t; Кх, Ку , Кг — коэффициенты турбулентной диффузии соответственно вдоль осей Ох, Оу, Ог; и — компонента средней скорости ветра вдоль оси Ох; ы — скорость осаждения частиц примеси вдоль оси Ог на подстилающую поверхность; г0 — коэффициент шероховатости подстилающей поверхности; р(х,у,г), /,

— соответственно фоновая концентрация, функция источника, скорость сухого осаждения этой примеси.

Соотношения (1)—(4) определяют математическую модель процесса рассеяния примеси в турбулентной атмосфере [3].

Цель данной работы — предложить метод определения коэффициента турбулентной диффузии Кг по экспериментально заданным значениям концентрации примеси q(t, х, у, г), мощности точечного источника непрерывного действия Q(t) и параметрам модели (1) — (4): и , ы , Кх , Ку .

Необходимость вычисления значений Кг по другим заданным значениям параметров математической модели (1) — (4) продиктована большими затруднениями, возникающими при экспериментальном определении его значений [3, 4].

(1)

q(0, х, у, г) = 0 ,

(2)

[Кг |г + wq}\2=20 = [У^}\г=го ,

(3)

q(t, х, у, г) ® 0, х2 + у2 + г2 г > 0,

(4)

2. Методика решения задачи определения вертикальной составляющей коэффициента турбулентной диффузии

Согласно [4] коэффициенты турбулентной диффузии Kx, Ky имеют

вид:

Kx = Ky = K0U , K0 = const, U = U(t, z).

Поэтому задача определения Kx и Ky сводится к задаче определения

U. Последняя - не вызывает на практике больших затруднений, поскольку современными техническими средствами легко определить изменения U от времени t и координаты z. Основная трудность заключается в нахождении коэффициента Kz(t,z) .

Пусть источник f в (1) является точечным с координатами (x0, y0, H), т.е. [3]

f (t, x, y, z) = Q(t)8(x - X0)d(y - ,0)8(z - H), где 8(s - s0) - дельта-функция Дирака, Q(t) - количество примеси,

выбрасываемой источником в момент времени t .

Согласно [4] коэффициент турбулентной диффузии Kz (t, z) возрастает в приземном слое атмосферы пропорционально высоте z :

Kz = K1(t)z , (5)

где K1(t), t e [0,T], - согласно поставленной задаче, неизвестная функция подлежащая определению.

Из (5) и (1) следует, что

K1(t) =

^+u ^ - w -A kx dq-A k, ^

dt dx dz dx dx dy Эу

дд + Э2£ . (6)

Эг Эг2

Таким образом для решения рассматриваемой обратной задачи достаточно вычислить

Эq Эq Эq Э2q Э^ и Э2q Э^ Эх ’ Эг ’ Эх2 ’ Эу2 Эг2

в заданных точках (х, у, г) в момент времени г и подставить эти значения в правую часть (6).

Согласно [5] задача нахождения производной п -го порядка г(г) функции и (г) (т.е. г(г) = и (п)(г)) сводится к решению (относительно z(t))

интегрального уравнения первого рода. В частности, для г(г) = имеем

Эг

уравнение

| 2(т)а?т = и (г) - и(0),

(7)

для г(г) =

Э 2и(г) Эг2

- уравнение

} (г - т) • г(т)ёт = и(г) - и(0) - г

Эи (г)

V Эг М=0

(8)

Предполагаем, что и(0)

Эи(0)

Эг

заданные величины.

Обозначим

Кх(,х,у,2) = МхЫ>, к(г,х,у,г) = ШЪЫ1, К, (,х,у,2)=ШЪЫ1 ,

х

Эг

э,

Кхх(ху»г) = Э 4 (Э Уу ), Куу(ху,г) = Э ^(Эх;уг), ягг(х,у,г)= Э уг) .

Эх2 Эу2 Эг2

Тогда (см. (7),(8)) для определения, например, яг (г, х, у, г) и ягг (г, х, у, г) будем иметь интегральные уравнения:

г

д(г, х, у, г) - д(г, х, у,0) = | (г, х, у,т)^ . (9)

д(г, х, у, г) - д(г, х, у,0) = | [(г - т) • Я22 (г, х, у, т)]оГ + г • Яг (г, х, у,0).

(10)

Соотношения (9) и (10) представляют собой интегральные уравнения первого рода относительно неизвестных функций К и ягг соответственно. Задача построения решения таких уравнений является некорректно поставленной [4].

0

0

0

При решении этой задачи перейдем от (9), (10) к их дискретным аналогам [3]:

q (г,х,у,г) -q (г,х,у,0) = £[Я2(г,х,у,^) • гк], (11)

к=1

q (г,х,у,г) - q (г,х,у,0) - г • Яг (г,х,у,0) = £[(гр - ^)Ягг(г,х,у,^) • Гк], (12)

к =1

г1,...,гр - точки деления интервала [0, г],

Гк =

, к = 1,

2

гк+1 - гк-1

2

к = 2,..,(р -1), (13)

гр гр-1 ,

----—, к = р.

2

Согласно (11), (12) по значениям q(t1,х,у,г),..., q(tN,х,у,г), заданным в точке (х,у,г) в различные моменты времени t1,..,tN е[0,s] с ошибками измерения соответственно п1 = у(^), п2 = п~(,2) , пN = v(tN) (п~(,) - случайный

процесс типа белого гауссова шума), требуется найти (восстановить) значения Яг (,ь х, у, гк),..., Яг ^, х, у, гк ) и Ягг (,ь х, у, гк),..., Ягг ^, х, у, гк)

соответственно, к = 1,...,р .

Введем в рассмотрение матрицу А = (Ак), все столбцы которой одинаковы, для решения уравнения (11) матрица А имеет вид:

А1к = Гк 5 к = 1,...5 р 5 1 = N 5

для решения уравнения (12) - вид:

А,к = (гр - гк) • Гк, к = 1,..., р, I = 1,...,N .

С учетом введенных выше обозначений и замечаний из (11) получим следующую систему линейных алгебраических уравнений:

хкк • Яг(*1,X,у,)]+VI = q(tl,х,у,г) - д(^,х,у,0), к=1

X • • • [2, X, у, гк)]+V2 = ]2, X,у, г) - q(t2, X, у,0), к=1

X • • • ^, X, у, гк)]+VN = ]($я, X, у, г) - ^(^^, х, у,0), к=1

(14)

из которой надо определить Яг ^, X, у, гк ) .

Из (12) получим соответствующую систему линейных алгебраических уравнений:

X кк • Ягг (^, X, у, гк)]+V! = q(tl, X, у, г) - q(tl, X, у,0) - г • Яг (^, X, у,0), к=1

X • • • ^2, X,у, Ек)]+V2 = q(t2, X, у, г) - q(t2, X, у,0) - г • Яг ^, X,у,0),

к=1

(15)

X • • • ^, X, у, Ек)]+VN = q(tN, X, у, г) - q(tN, X, у,0) - г • Яг (^, X,у,0),

к=1

из которой надо определить ЯгЕ X, у, гк ) .

Систему (14), (15) представим в матричном виде:

ЛЯг + V = 'д , (16)

ЛЯее +у = 'д , (17)

где , X, у, Ек) = q(tl, X, у, Ек) - q(tl, X,у,0),

q(t1, X,у, гк) = q(tl, X,у, ^) - q(t1, X, у,0) - г • Яг (^, X, у,0), I = (1,2,...,N), к = (1,2,...,р).

Для подавления влияния значений белого шума ) на значения

Яг (к X,у,гр) и Ягг (tk, X,у, гр) , к = 1,...^, можно использовать многошаговый

(многократный) фильтр Калмана-Бьюси [6].

Для этого задаем начальные приближения для решения яЕ0) = яг (0, X, у, г) и матрицы ковариаций ошибок решения Р(0). Для их выбора удобно использовать метод регуляризации Тихонова [4], согласно которому

я.

(0)

= (аЕ + ЛтЛ) 1 ЛTq , Р(0) = 82 (аЕ + лтл)~1,

(18)

где Е - единичная матрица, а > 0 - параметр регуляризации, играющий роль неопределенного множителя Лагранжа, 5 > 0 - верхняя оценка значения погрешности правой части (16).

Последующие приближения яЕ1) решения яЕ системы (14) могут быть найдены по следующей итерационной схеме [6]:

Зададим начальные приближения для решения 40) = ягг (0, X, у, Е) и матрицы ковариаций ошибок решения р(0). Для их выбора удобно

использовать соотношения (18), подставив в них я^ вместо дЕ°).

Последующие приближения я£) решения ягг системы (15) могут быть найдены по итерационной схеме (19)-(20) путем замены я£) на яЕ1).

На практике можно столкнуться с ситуацией, когда обратные матрицы в соотношениях (18)-(20) найти (определить) невозможно (рассматриваемые матрицы могут быть вырожденными). В этом случае вместо обратных матриц следует использовать в (18)-20) псевдообратные, воспользовавшись методом Гревилля построения псевдообратной матрицы [7].

Соотношения (18)-(20) позволяют найти значения величины я^ -оценку ЯЕЕ с заданной погрешностью е > 0. Способ нахождения оценки яЕЬ) для яЕ также подробно описан. Аналогично определяются я^, яУу, я(ь), яХУ) соответственно для яхх, яуу , Я{ , ях .

Подставляя найденные оценки в (7), получим наилучшую в среднем квадратическом смысле оценку к^к) значения К^к):

V

/

(19)

р{1) =ИрА-1) ) + ((1)) х(1)А(1) 1 , ж(1) = М[у(1) (~(г) )Г], I = 1,2,..., у . (20)

V

/

К«к) =

+

(Ь)

■ку • Куу

+ -

дКУ. К(Ь). ду у

К ■ я

(Ь)

К

дг

■К

(Ь)

б('к)

ЯЬ) + Я(Ь)

(20)

3. Пример

Для проверки качества работы алгоритма по указанной методике, воспользуемся экспериментальными данными, взятыми из отчетов Центра лабораторного анализа и технических измерений по Южному Федеральному округу (ЦЛАТИ по ЮФО) и содержащими информацию о выбросах в атмосферу диоксида азота. Согласно этим данным:

б = у+1 (кг/с), Н = 20 м, и = 0,51п г (м/с), Кх = Ку = К0и м2/с, К0 = 0,25 м, г0 = 0 с,

^ = 0,01 (м/с).С помощью (20) найдены наилучшие в среднем квадратическом смысле оценки значения вертикальные составляющие коэффициента турбулентной диффузии на промежутке времени г е [0;55] (вычисления проводились в пакете прикладных программ ММЬаЬ). Графическая визуализация результатов проведенных расчетов приведена на рисунке 1.

Кг О

8

6

4

2

О

О 10 20 30 40 50 г

Рисунок 1 - Графическое изображение совпадения значений

экспериментальной и расчетной вертикальной составляющей

коэффициента турбулентной диффузии.

СПИСОК ЛИТЕРАТУРЫ

1. Алоян А.Е. Моделирование динамики и кинетики газовых примесей и аэрозолей в атмосфере. - М.: Наука, 2008. - 415 с.

2. Семенчин Е.А., Кармазин В.Н., Калина Н.Н. О разрешимости некоторых обратных задач для уравнения атмосферной диффузии. Экологический вестник научных центров Черноморского экологического сотрудничества, №4, 2005. - С. 47-51

3. Семенчин Е.А. Аналитические решения краевых задач в математической модели атмосферной диффузии. Ставрополь: СКИУУ, 1993. - 141 с.

4. Берлянд М.Е. Современные проблемы атмосферной диффузии и загрязнения атмосферы. - Л.: Гидрометеоиздат, 1975. - 448 с.

5. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М.: Наука. Главная редакция физико-математической литературы, 1979. - 142 с.

6. Сизиков В.С. Устойчивые методы обработки результатов измерений. Учебное пособие - СПб: Изд-во «СпецЛит», 1999. - 240 с.

7. Гантмахер Ф.Р., Теория матриц. - Москва: изд-во "Физматлит", 2004. - 576 с.