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

ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ

«НАУКА. ИННОВАЦИИ. ТЕХНОЛОГИИ», №3, 2013

УДК 519.5 Наац В. И. [Naats V. I.],

Гаршина Т. В. [Garshina T. V.]

ВЫЧИСЛИТЕЛЬНАЯ МОДЕЛЬ НЕСТАЦИОНАРНОГО УРАВНЕНИЯ ПЕРЕНОСА ПРИМЕСЕЙ НА ОСНОВЕ МЕТОДА ВЗВЕШЕННОЙ НЕВЯЗКИ И ОПЕРАТОРОВ ОБОБЩЕННОГО ДИФФЕРЕНЦИРОВАНИЯ ФУНКЦИЙ

The calculation model for the unsteadily-state mass transfer equation, which is based on weigh discrepancy method and the operators of the generalized differentiation of the functions

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

Ключевые слова: нестационарное уравнение переноса примесей в атмосфере, рекурсивный вычислительный алгоритм, сингулярные интегралы, численные исследования.

The non-stationary equation of transfer of impurity in the atmosphere is considered, which involves the use of empirical data. Construction of the corresponding recursive algorithm is performed on the basis of weigh discrepancy method and the operators of the generalized differentiation of the functions.

Key words: non-stationary equation transfer of impurity in the atmosphere, a recursive computing algorithm, singular integrals, numerical researches.

Представленная работа продолжает исследования авторов по разработке численных методов и алгоритмов для дифференциальных уравнений в частных производных в ситуациях, когда исходные данные являются приближенно заданными функциями [1, 2, 3]. Подобным примером является задача прогноза экологических ситуаций в пограничном слое атмосферы, в пределах которого распространение и рассеяние загрязняющих веществ описывается так называемым уравнением перено-

са. В этом случае исходные данные представляются приближенными значениями полей скорости ветра V(P, t) и коэффициента турбулентного обмена K(P, t). Если значения поля V(P, t) в точке P(x, y, z) в определенные моменты времени t могут быть оценены прямыми измерениями в пределах контролируемого региона, то значения функции K(P, t) вводятся в уравнение переноса на основе статистических данных либо с помощью полуэмпирических формул. Постановка задачи, о которой пойдет речь ниже, формулируется с учетом сделанных выше замечаний.

Допустим, что концентрация загрязняющих веществ в пограничном слое атмосферы описывается распределением u(x, t), если рассматривать одномерный вариант задачи (скажем, перенос субстанции в направлении координатной оси Ox). Тогда исходное функциональное уравнение, связывающее поле концентрации u(x, t) с характеристиками реальной среды V(x, t) и K(x, t) представляется в виде следующего дифференциального уравнения параболического типа:

— (x, t) = — (K(x, t)u') - — (V(x, t)u) + S(x, t), (1)

dt dx x dx ' w

где (x, t) € О = Ox xfit = [a,b] x [tg,T], S(x, t) - функция источ-

ника загрязняющих веществ.

Для определенности решений уравнения (1) задаются так называемые начальные и граничные условия u(x, t0) = u0 (x), u(a, t) = u1 (t), u(b, t) = u2 (t). Значения полей V(x, t) и K(x, t) в области W образуют множество исходных данных рассматриваемой задачи, которые далее лучше обозначать как B{V(x, t), K(x, t)} .

В соответствии с уравнением (1) от функций V(x, t) и K(x, t) требуется их дифференцируемость по переменной x в пределах Wx, т. е. существование непрерывных частных производных Vx(x, t) и K'x(x, t). Это означает, что элементы множества B должны принадлежать классу функций C1 (W). Если же данные V(x, t) и K(x, t) формируются на основе эмпирических наблюдений, то их следует считать приближенными и не дифференцируемыми в обычном смысле функциями. В соответствии с этим вводим обозначения B\V(x,t),K(x,t)} и считаем B с С^ (п.), т. е. элементы множества B в лучшем случае суммируемые функции на W. В дальнейшем в целях большей ясности и определенности полагаем, что указанные

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

IIIX - ¡а(х)11т)< а III , (2)

где s - достаточно малое число (s > 0).

Поскольку функция f(x) обычно неизвестна и подлежит определению, то в оценках типа (2) можно использовать приближение II f(x)\\h)f \\h). Отличие функций f(х) и f(x) состоит в том, что для f (x) в принципе применим оператор обычного дифференцирования D = djdx по переменной x , в то время как для «дифференцирования» fs (x) требуется построение специального оператора Ds в той ли иной мере аналогичного исходному оператору D , определенному на множестве CДП) . В этом контексте оператор Ds называется оператором обобщенного дифференцирования, и областью его определения являются более широкие функциональные классы, нежели CДЦ,) . Построение оператора Ds может быть осуществлено различными методами, выбор которых согласуется с особенностями решаемых прикладных задач, но в любом случае требуется, чтобы при s ® 0 оператор Ds в том или ином смысле был близок к оператору L. Напомним, что последний определяется через предельное отношение (А/(x) / Ах) при Ax ^ 0 .

Если f (x) не дифференцируется в некоторой точке x, то применение конечно-разностных аппроксимаций для оператора D мало обосновано и практически ведет к расходимости вычислительных процессов. В пределах данной работы понятие обобщенного дифференцирования связывается с построением решающего алгоритма для модели (1) в рамках вариационного подхода, в том смысле как он понимается в работах [1, 4, 5, 6]. Дадим соответствующие пояснения по существу вопроса.

В основе подхода лежит предположение о возможности представления исходного решения u(x, t) в следующем виде

m m+1

u{X,t) ^ Um(x,t) - Zck(t)vk(x) + C0(t)v0(x) + Cm+i(t)Vm+i(x) - ^Ck(t)vk(x), (3)

k-1 k-0

где {vk (x )} - некоторая система базисных функций.

Подробнее о выборе базиса для моделей типа (1) можно найти в работе авторов [1]. Обычно на функции базиса в краевых задачах

налагают Условия: Vk {Х ) I =а = Vk {Х ) \х = = 0 , V0{X ) \х =а = Vm+1(x ) I= = 1 . В результате в выражении (3) C0(t) = u1(t), Cm+1(t) = u2(t) и, значит, в разложении (3) искомой функции неизвестными остаются коэффициенты C1(t) ,...., Cm(t) (тоже вектор-функция C(t) ). Остается построить систему функциональных уравнений для определения указанных коэффициентов. С этой целью исходное уравнение (1) перепишем в виде

д

U(x, t) =--[V(x, t)u - K(x, t)u'x ] + S(x, t) . (4)

dx

Введем обозначение g(x, t) = u1 (t)v0 (x) + v2 (t)vm+1 (x) и подставим представление искомой функции u(x, t) (3) в (4). После соответствующих преобразований придем к следующему дифференциальному уравнению относительно компонент вектора - коэффициентов C :

m m О

£ Ck (t )vk (x ) + g(x, t) = Ck (t) д- [V (x, t )vk (x ) - K (x, t)v'x (x )] -

k=1 k=1 dx

д

д- [V (x, t)g(x, t) - K (x, t )g'x (x, t )] + S(x, t) . (5)

Уравнение (5) можно рассматривать как некий приближенный аналог исходного уравнения (1) в предположении, что представление (3) для искомой функции u(x,t) действительно имеет место. Здесь имеется ввиду, что um (x, t) ® u(x, t) при m ^ ж для каждой пары точек (x, t) € О . Вместе с тем следует заметить, что подобная эквивалентность (1) и (5) более чем условна, ибо из близости um(x,t) и u(x,t) никаким образом не следует близость u'mx(x, t) к u'x(x, t), не говоря уже о вторых частных производных. Положение усугубляется, если функции V(x,t) и K(x,t), входящие в выражение (5), заданы приближенно и строго говоря не дифференцируемы.

К этому вопросу вернемся позже, а пока рассмотрим задачу построения системы m уравнений для определения m неизвестных функций Ck(t), (k = l,..m), исходя из функционального уравнения (5). Ясно, что предварительно следует надлежащим образом «свернуть» уравнение (5)

по переменной x. Подобная операция в рамках вариационного подхода осуществляется интегрированием невязки p(x, t, C) для уравнения (5) по переменной x, требуя при этом выполнения условия:

^wt(x)p(x, t,C)dx = 0 , l = l,..m , (6)

nx

где {w (x)} - система весовых функций.

Вариантом выбора {wl (x)} является случай, когда весовые функции w (x) локализованы в пределах частичных (элементарных) интервалов Д1 (x) = xt _ xt—1, покрывающих область = [a,b]. В этом случае пару (А; (x), w{ (x)) называют конечным элементом, а подход в целом -методом конечных элементов [7]. Применяя операцию (6) к уравнению (5), приходим к следующей системе линейных дифференцированных уравнений относительно компонент Ck (t), (k = 1,. .m)

m m

ECk = -ECk (t\(t) + ъ, (t), (7)

k=1 k=1

где обозначено:

b

aik = f wi (xH (x)dx , k,1 = 1m , (8)

bk = f W X faJ (X't)dX (9)

а

Jk(x, t) = J(x, t, vk) = V(x, t)vk(x) - K(x, t)v'(x), (10)

b

h(t) = J w(x){g(x, t) + J(x, t, g) - S(x, t)}dx . (11)

a

Представленная выше система (7)-(11) может быть записана в матричном виде

AC = -B(t)C + h(t) . (12)

Численные методы решения подобных систем, а так же соответствующие программные средства, применяемые для решения задач переноса субстанции, описаны в работах авторов [1, 5, 6]. Здесь лишь заметим, что одним из простейших методов является построение рекурсивной схемы вида

с( з +1) = т(з )с(з) + тб(])Л(]),

где оператор шага = (I + тА^В^))'1^ - тА^В^)) и

оператор источника §(]) = (I + ;тА~1Б(]))~1. В этих выражениях параметр т = Ь.+1 — Ь. - интервал дискретизации искомой функции и(х, Ь) по переменной Ь.

Теперь вновь обратимся к проблеме дифференцируемос-ти функций, которые так или иначе задействованы в вычислительной схеме (7)-(11). Ясно, что если исходные данные V(х, Ь) и К(х, Ь) принадлежат множеству Ва, то вычисление Л'кх(х, в выражении (9) затруднительно. Положение спасает то обстоятельство, что указанная производная входит под знак интеграла (как следствие операции усреднения [6]). В связи с этим к интегралу (9) применим формулу интегрирования по частям и получим следующее выражение:

ь ь

J (х(х, г)йх = ■ш,/к (х, г)йх + [ш1 (х)Лк (х, . (12)

а а

Ясно, что необходимо выполнение предположения о дифференцируемости весовых функций в той же мере, как выше предполагалась дифференцируемость функций исходного базиса Ук (х) (см. [10]). Интегральное равенство (12) делает вычисление интеграла (9) вполне определенным для любой пары (V , Ка) е Ва , ибо требует от соответствующих функций лишь их интегрируемости. Ясно, что при а ® 0, когда (V , Ка) ® (V, К) в каждой точке (х, Ь) € О, формула интегрирования по частям (12) обратима и, значит, оператор , стоящий за равенством (12), эквивалентен оператору обычного дифференцирования О. Последним объясняется смысл термина «обобщенное дифференцирование».

Располагая описанным выше вычислительным алгоритмом, можно непосредственно обратиться к проблеме надлежащего выбора весовых функций, сообразуясь с теми требованиями, которые оговаривались выше в процессе построения (7)-(11). Эти требования могут быть сформулированы следующим образом: а) (х) > 0 для всех х € (I = 0,1..т + 1); б) /'Ш1(х№ = 1 для У1 = 1,т; в) система {т1 (х)} (I = 1,..т) может быть ассоциирована с некоторым множеством узловых точек {х(} (I = 0,..(т + 1)), расположенных на отрезке = [а,Ь] таким образом, что (х) = Цж, х 1); г) функция -ш(х, х') такова, что т(х, х') = и)(х — х'), (х, х') € хОх,, и при этом Цх, х') ® 0 при | х — х' ^ ; д) производная (х, х') суммируема на Нх по переменным х и х'.

Следует заметить, что подобные функции известны в прикладном анализе, и в частности используются в конструктивной теории функций и теории приближений [8]. Речь идет прежде всего о ядрах так называемых сингулярных интегралов функций. Напомним, что сингулярный интеграл для суммируемой функции, скажем, f (х), определенной на отрезке [а, Ь] записывается следующим образом:

/ К (х, х')/(х ')<Ъ' = (К/) (х) = / (х), 11т / (х) = /(х),

(13)

где

Кп (х, х') - ядро интеграла (13) порядка п (п ^ х>).

В теории сингулярных интегралов [9] последовательность ядер {Кп (х.х')} удовлетворяет всем перечисленным выше требованиям на весовые функции в интеграле (6). Определяющее условие (13) для сингулярных интегралов (Кп/)(х) распространяется на производную f '(х), если последняя является интегрируемой функцией. Имеем

/ К (х, х')/ '(х >х' = (К о Б/)(х) = (Б/)п (х), Ит(БЦ (х) = (Б/ )(х).

(14)

Если применить к интегралу (14) формулу интегрирования по частям, то придем к следующему равенству

f Kn (x, x')f(x')dx' = -f K[x (x, x ')f(x ')dx'+ [Kn (x, x y(x)l

При n [Kn (x, x')/(x ')]XЦ ^ 0, и учитывая, что

K,x' = найДем

б б f Kn (x>x')/(x l)dx' - / Klxx (x>x')/(x Vх'

a

a

/

a

В силу предельного соотношения в (14), можно утверж-

b

дать, что интеграл / К,(x,x')/(x ])dx' является приближением к f '(x) для дифференцируемых функций. Именно о последних функциях и идет речь в численных методах для дифференциальных уравнений.

Ясно, что интегральный генератор с ядром (DxKn )(x, x') = K'nx (x, x') определяет отображение f (x) ® f '(x) в каждой точке x, в которой f (x) дифференцируема.

Поскольку подобное отображение формально определено для всякой суммируемой функции f (x), то этот оператор можно считать оператором обобщенного дифференцирования. Для недифференцируемых функций f(x) значения интеграла (Dx о Knf) (x) в точке x характеризуют меру гладкости f(x) в окрестности этой точки. С учетом вышеизложенного формальная операция, связанная с применением формулы интегрирования по частям к функции Jk(x, t) по переменной x (см. [12]), приобретает заметно больший содержательный смысл при условии w{ (x) = w(x, xt) = Kn (x, xt). Более полное изложение теории интегральных представлений искомых функций типа (13) приведено в работах [1, 2, 3].

В заключение остановимся кратко на выборе последовательностей {Kn(x,x')} в прикладных задачах. Ясно, что решение подобных вопросов целесообразно согласовывать (в пределах возможного) с особенностями исходных функциональных уравнений параболического типа, каковым является модель (1). Заметим, что во всех приложениях, когда решается уравнение типа U = auxx + s(x, t), то важную роль играет функция

E(x, t) = expI- — I, (a > 0), (15)

2Vpat I 4atI

известная как фундаментальное решение однородного уравнения указанного типа. Функция Е(х, Ь) позволяет решение неоднородного уравнения представить в виде следующего интеграла:

и(х, Ь) = J J Е(х - хь - ь '>(хь х Ш'. (16)

Это интегральное представление во многом предопределяет аналитические свойства функций, образующих множество решений для неоднородных уравнений параболического типа. Функция Е(х, Ь) удовлетворяет условию нормирования, а именно:

Г —^=ехР\- —— г' i [ j

для всех t > 0 .

Это позволяет на основе Е(х, Ь) в данном классе задач сконструировать так называемое распределение единицы. Действительно, выражение (15) нетрудно привести к виду

Wя (х, X ',<) = -П= ехр П2{Х- Х 1)21, (17)

где обозначено d(t) = 44м и п = 1,2...

Нетрудно показать, что функции последовательности Шп (х, х ) удовлетворяют ограничениям, налагаемым условиями (а)-(д) на функции т(х, х'). В дальнейшем удобно ввести параметр т = уп и формулу (17) писать в виде

W (x, x ',d) = — exp \nd ■ т

\2

x — x'

d т

, (0 < t < 1).

Условие нормировки для (16) выглядит следующим образом:

b

lim J Wt (x, xd)dx1 = 1, "(x,d),

а

где x E (a,b), d > 0 . d > 0

Теперь в методе взвешенной невязки принимается соотношение:

W (x):= Wt {x{, x,d), l = 1, m, (18)

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

С учетом (18) расчетные формулы (8), (9) и (11) решающего алгоритма примут следующий вид

b

а(Т)(t) = f W (x,, x, t)vk (x)dx,

a b

W = f K,x X , x, t)Jk (x, t)dx + t) ,

a b

h!¡r)(t) = f W(x¡, x, t){g(x, t)Jk(x, t, g) - S(x, t)}dx ,

a

Итоговая система уравнений теперь имеет вид

C = -A-lBC + A-lh , (19)

где операторы A— и B являются операторозначимыми фун-

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

В конкретных приложениях изложенной здесь теории на основе вычислительного эксперимента подбирается функциональная зависимость т(а) такая, что r(a) ® 0 при s ® 0 . Решением системы (19) является вектор C(т), определяющий согласно (3) последовательность приближенных решений |и(т'(x, t)} по параметру т . Последующая оптимизация вычислительных алгоритмов по параметрам задачи таким, как

параметры т, а , размерность задачи т, интервал временной дискретизации Т и других, может быть осуществлена на основе принципа мини-макса [10].

ЛИТЕРАТУРА 1.

8. 9.

10.

Наац В. И., Наац И. Э. Математические модели и численные методы в задачах экологического мониторинга атмосферы: монография. М.: Физматлит, 2010.

Рыскаленко Р. А., Черкасова И. В. Интегральные представления функций в численных методах решения нестационарных задач переноса // Вестник СКФУ. 2013. № 1 (34). Ставрополь: Изд-во СКФУ, 2013. С. 30-34.

Рыскаленко Р. А., Чемеригина М. С. Операторы обобщенного дифференцирования в численных методах решения нелинейного уравнения переноса с приближенными данными // Вестник СКФУ. 2013. № 1 (34). Ставрополь: Изд-во СКФУ, 2013. С. 35-38.

Ректорис К. Вариационные методы в математической физике и технике / пер. с англ. М., 1985.

НаацВ.И., ТравкинаТ.В. Разработкарекурсивноговычислительного алгоритма для оценки атмосферной турбулентности на основе уравнения переноса и результаты вычислений // Вестник СКФУ. 2013. № 1 (34). Ставрополь: Изд-во СКФУ, 2013. С. 9-14.

Наац В. И. Вычислительная модель нестационарного уравнения переноса примеси на основе метода взвешенной невязки // Известия вузов Сев.-Кав. региона. Естеств. науки. Прил. 5'04. Ростов-н/Д, 2004. С. 3-15.

Зенкевич О., Морган К. Конечные элементы и аппроксимация / пер. с англ. М., 1986.

Натансон И. П. Конструктивная теория функций. М.: Физматлит, 1949.

Данфорд Н., Шварц Дж. Т. Линейные операторы: общая теория. М.: Изд-во иностранной литературы, 1962.

Сухарев А. Г. Минимаксные алгоритмы в задачах численного анализа. М., 1989.

ОБ АВТОРАХ Наац Виктория Игоревна, ФГАОУ ВПО «Северо-Кавказский федеральный университет», доктор физико-математических наук, профессор кафедры математического анализа, тел. (8652) 35-2110,

E-mail: VINaac@yandex.ru.

Гаршина Татьяна Васильевна, ФГАОУ ВПО «Северо-Кавказский федеральный университет», соискатель кафедры математического анализа,

E-mail: travkinatv@mail.ru.

Naats Viktoria I., North-Caucasian Federal University, doctor of physical and mathematical sciences, Professor of the Department of Mathematical analysis.

Garshina Tatiana V., North-Caucasian Federal University, competitor of the Department of Mathematical analysis.