Аналитическое и численное построение решений системы уравнений газовой динамики, имеющих спиральный характер Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 16, № 1,2011

Аналитическое и численное построение решений системы уравнений газовой динамики,

о >к

имеющих спиральный характер*

С. П. Баутин, А. В. Рощупкин Уpальский госудаpственный университет путей сообщения,

Екатеринбург, Россия e-mail: [email protected], [email protected]

Рассматриваются такие решения системы уравнений газовой динамики, которые передают плоские изoэнтропические течения политропного газа при учете силы Корио-лиса, вызванные заданным стоком газа на окружности ненулевого радиуса. Доказано, что в начальные моменты времени в исходном однородном покоящемся газе одновременно со стоком возникает закрутка газа: в положительном направлении в Северном полушарии и в отрицательном в Южном. С помощью интегрирования системы обыкновенных дифференциальных уравнений получено стационарное решение системы уравнений газовой динамики, передающее течение с достаточно большой закруткой газа в окрестности окружности стока и имеющее нулевую закрутку на окружности заданного притока газа. Численным методом характеристик построено нестационарное решение системы уравнений газовой динамики, описывающее плоское спиральное течение с заданным стоком.

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

1. Система уравнений газовой динамики при наличии силы Кориолиса

Система уравнений газовой динамики (СУГД) для изoэнтропических плоских течений политропного газа при учете силы Кориолиса имеет следующий вид [1]:

Ct + UCr +

(7-1)

о,

< Ut + UUr — JT +

(7-1)

ccr = av,

(1

Uv

„ Vt + uvr + -Црг- = — au.

Здесь c = p(l-1^/2 — скорость звука; y = const > 1 — показатель политропы газа в уравнении состояния p = pY/y; в плоскости переменных x, y введена полярная система координат (r,p) и предполагается, что д/дp = 0; u, v — соответственно радиальная и окружная составляющие вектора скорости газа; a = 2П sin ф — параметр Кориолиса; П — модуль угловой скорости вращения Земли; ф — широта точки O на поверхности Земли,

"Работа выполнена при финансовой поддержке РФФИ (грант № 08-01-00052).

2

в которой находится начало координатной плоскости xOy, касающейся поверхности Земли в точке O. Если точка O лежит в Северном полушарии, то 0 < ф < п/2, в Южном — —п/2 < ф < 0, на экваторе — ф = 0.

В системе (1) стандартным образом введены безразмерные переменные. Если в качестве масштабов скорости и расстояния при введении безразмерных переменных взяты

соответственно с0о = -103м/с, г0о = Ю3 м, то безразмерное значение константы О =

3

0.000218, а ¿оо — масштабное значение времени — равно трем секундам. Если, не меняя масштаба скорости, увеличить масштаб расстояния, например в сто раз, то значения констант О и t00 также возрастут в сто раз.

2. Решение в начальные моменты времени задачи о плавном стоке газа

Рассматривается следующая задача о плавном стоке газа.

Пусть задана некоторая аналитическая функция u0(t) такая, что

ио(0) = 0, u0(0) < 0. (2)

Для системы (1) ставятся условия

c(t,r)\c+ = 1, u(t,r)\c+ = 0, v(t,r)\o+ =0, u(t,r)|r=ro = U0(t). (3)

Решение задачи (1), (3) при условиях (2) описывает плавный сток на окружности r = r0 > 0 первоначально однородного и покоящегося при r > r0 газа.

Схема течения в задаче о плавном стоке представлена на рис. 1, где цифрой 0 обозначена область покоящегося газа, цифрой 1 — область определения решения задачи (1), (3).

Первые три условия из (3) обеспечивают непрерывное примыкание решения задачи (1), (3) через звуковую характеристику

C + : r = r0 + t

к однородному покоящемуся газу. Функция u0(t) из четвертого условия в (3) обеспечивает единственность решения поставленной задачи и задает на окружности r = r0 скорость стока газа, непрерывно уменьшающуюся с ростом времени от нулевого значения в момент времени t = 0.

На рис. 2, a приведено качественное поведение функции u0(t) в случае плавного стока. Теорема. Задача (1), (3) при условиях (2) имеет как в некоторой окрестности точки (t = 0,r = r0), так и в некоторой окрестности C+-характеристики во все моменты времени единственное аналитическое решение.

Рис. 1. Схема течения в задаче о плавном стоке газа

Рис. 2. Качественное поведение функции и0(Ь) и газодинамических параметров и(Ь, т)|^=С0^ (б), с(Ь, т) |^=С0пэ! (в) и у(Ь, т) |^=С0пэ! (г) в случае плавного стока в моменты времени Ь = Ь\}2 при Ь\ < Ь2 (Северное полушарие); 1,2— кривые, соответствующие моментам времени Ь\,Ь2

Данная теорема справедлива, так как задача (1) — (3) является частным случаем характеристической задачи Коши (ХЗК) стандартного вида [2]. Сведение задачи (1) — (3) к нужному виду фактически повторяет соответствующие выкладки из работ [3—5] и поэтому здесь не приводится.

Решение ХЗК (1), (3) представляется в виде сходящегося ряда

к=0 к-

д к и

дгк

(4)

с+

где

И =

/ с N

и

\ V )

есть вектор искомых функций, а И0 (Ь) — постоянный вектор начальных данных на характеристике С +, составленный из первых трех равенств в (3).

Построение коэффициентов ряда (4) при к > 1 проводится стандартным образом [2] так же, как это делается при построении решения задачи о плавном движении поршня [3—5]. В частности, для определения ик (Ь), к > 1, следует решать обыкновенные дифференциальные уравнения с начальными данными, определяемыми в момент времени Ь = 0,с помощью четвертого соотношения из (3). После этого ск (Ь), ук (Ь) находятся в явном виде из соответствующих алгебраических уравнений. Детальный анализ структуры коэффициентов ряда (4) повторяет соответствующие построения из [3—5] и приводит к выводу о неограниченности области сходимости ряда по переменной Ь: при любом значении Ь > 0 имеется ненулевой радиус сходимости ряда (4) в окрестности С+ -характеристики, который, естественно, стремится к нулю при Ь —> +оо.

Первые коэффициенты ряда (4) имеют вид

1

ui(t) =-F-, , -^-- > О,

VtTñ¡

(7+l)v/^(J1 +1-1]- 1

Го ) ^/ñ¡u'0(Q)_

(y -1)

с\ = —-—щ > 0, v\{t) = 0, V2 = au\. (5)

Из последнего равенства в формулах (5) следует, что если течение рассматривается в Северном полушарии, где a = 2Ü sinф > 0, то v2(t) > 0, а если в Южном, то v2(t) < 0.

Поскольку доказано, что решение задачи (1) — (3) задается аналитическими функциями, то для коэффициентов ряда (4) справедливы стандартные оценки [6, 7]. В частности, и для коэффициентов ряда, задающего функцию v(t,r), существует неотрицательное число M такое, что при всех k > 0 в некоторой окрестности точки (t = 0,r = r0) справедливо неравенство \vk(t)(r — r0 — t)\k < Mkk!. Из этого неравенства следует, что в соответствующей окрестности рассматриваемой точки начальный отрезок ряда для v(t, r) (напомним, что v0 = vi = 0) задает главную часть функции v(t,r):

те

v2(t)(r — ro — t)2/2\ \vk(t)(r — ro — t)k/k!

k=3

поэтому под действием силы Кориолиса в течении, образующемся при стоке газа на окружности r = r0, при t > 0 в некоторой окрестности звуковой C + -характеристики начинается закрутка газа в положительном или отрицательном направлениях для Северного или Южного полушария соответственно.

На рис. 2, б—г для случая Северного полушария приведено качественное поведение газодинамических параметров в задаче о плавном стоке в разные моменты времени t = ti}2: ti < t2.

3. Стационарное плоское спиральное течение со стоком

У системы (1) при заданных значениях

r rin const > c\r=rin -L) u\r=rin uin) v\r=rin 0) где константа uin удовлетворяет неравенствам

— 1 < Uin < 0)

существует единственное стационарное (т. е. при d/dt = 0) решение

c = co(r), u = uo(r), v = vo(r),

определенное на промежутке

0 < ro < r < rin.

При этом функция vo(r) находится из третьего уравнения системы (1) при d/dt = 0 в явном виде

v°(r) = у - \т,сх = a-r¡n. (6)

По заданной функции п°(т) функция с°(т) определяется из первого уравнения системы (1) при д/дЬ = 0 равенством

с°(г)

Со

тп°(т)

(7-1)/2

Со

тгп Пг/

(7)

а функция п°(т) < 0 задается в неявном виде с помощью уравнения

2 / С \ (^-1) С2 П2

0,

(8)

являющегося следствием второго уравнения системы (1) при д/дЬ = 0 и полученных соотношений (6), (7). Здесь

Сз =

2

4- II --г

I I 4 ' 1.Т,

(7 _ 1) 1 ' 4 гп•

Условие существования неявно заданной функции

дЕ

дп°

= 0

может быть представлено с помощью конкретного вида частной производной

ди°

С2

тп°

(7-2)

С2

т(п°)2

2 и° = — п°

М2 " ( ^

тп°

(7-1)"

(9)

Выписать в случае произвольного значения константы ^ > 1 функцию п° = п°(т) в явном виде из уравнения (8) вряд ли возможно. Только в частном случае 7 = 3 (8) переходит в биквадратное уравнение относительно п°.

Функции с = с°(т), п = п°(т) можно также определить при решении задачи Коши

(с°)' = _

(7-1)

(и0)2 + $ - £ г2

т[(п°)2 _ (с°)2)

(п°)' = п°

(го\2 , С\ а2 2 Кс ) + -р: - 1"г

(10)

т[(п°)2 _ (с°)2]

.. с \т=Г1п 1, п \т=Г1п пгп-

Естественно, что определение функций с°(т), п°(т) с помощью равенства (7) и разрешения уравнения (8) эквивалентно решению задачи Коши (10), поскольку из (9) с учетом (7) следует

дЕ

2

дп° п°

2

(п°)2 _ (с°)

Решение задачи (10) (как и (6) — (8)) описывает следующее плоское стационарное спиральное течение газа при учете силы Кориолиса.

На окружности т = тгп осуществляется приток газа извне (из области с т > тгп), поскольку на ней задано п°(тгп) = пгп < 0 — отрицательное значение радиальной скорости газа. На этой же окружности заданы значения двух других газодинамических параметров:

о

Рис. 3. Графики функций с°(г) (а), и°(г) (б) и V°(г) (в), являющиеся решением одной задачи Коши (10); линии тока течения и окружности г = г0, г = гп, на которых происходит сток и приток газа (г)

с° (тп) = 1 — единичной скорости звука и V° (тп) =0 — нулевого значения окружной компоненты скорости.

Неравенства —1 < ип < 0 отражают дозвуковой характер течения в окрестности окружности притока т = тп и обеспечивают разрешимость как уравнения (8) относительно и° (т), так и задачи (10).

Плоское стационарное течение (6)—(8), т. е. решение задачи (10), определено в некотором кольце 0 < т0 < т < тп. Внутри этого кольца функции и°(т), V°(т) строго монотонны: и°(т) растет от и°(т0) до и°(тп), V°(т) убывает от V°(т0) > 0 до V°(тп) = 0. Во всем течении выполняется неравенство —и°(т) = с°(т).

Поскольку при 0 < т <т0 течение газа не строится, то на окружности т = т0 предполагается сток газа.

На рис. 3, а—в представлены графики функций с°(т), и°(т), V°(т), являющихся решением одной задачи Коши (10), на рис. 3, г приведены четыре линии тока течения и окружности т = т0, т = тп, на которых происходит соответственно сток и приток газа.

Таким образом, стационарные спиральные течения, имеющие заданный радиальный приток на окружности т = тп (т. е. и°\г=Пп = ип < 0, V°\г=пп = 0) и не имеющие особенностей на некотором отрезке 0 < т0 < т < тп, обладают следующими свойствами:

1) можно полагать, что на окружности т = т0 имеется сток газа;

2) на отрезке [т0, тп] скорость звука не является постоянной и при т — т0 + 0 монотонно убывает;

3) при т0 < т < тп радиальная скорость имеет отрицательные значения и при уменьшении т ее модуль растет;

4) окружная скорость определена (см. (6)) при всех положительных значениях т, равна нулю при т = тп и при т — +0 для случаев соответственно Северного или Южного полушария неограниченно растет или убывает;

5) при т — т0 + 0 численные значения модуля окружной скорости на порядок больше модуля радиальной скорости.

Поскольку при г ^ г0 + 0 величины \п°(г)| растут, а значения с°(г) убывают, то можно предположить, что для конкретного набора констант 7, а, гп, щп < 0 найдется ненулевое значение г = г* > 0, при котором [п°(г*)]2 = [с°(г*)]2, поэтому данное стационарное спиральное течение с притоком на окружности г = гп нельзя продолжить при г < г* и, следовательно, сток газа заведомо надо будет задавать на окружности с радиусом г0, где 0 <г* <го < гы.

4. Численное построение нестационарного спирального течения со стоком

Доказанная выше теорема обеспечивает существование нестационарного течения со стоком в некоторой окрестности точки (£ = 0,г = г0), а также при £ ^ в некоторой окрестности звуковой С + -характеристики, т. е. в моменты времени, близкие к моменту начала плавного стока из однородного покоящегося газа, сопровождающегося возникновением закрутки газа, а также в окрестности С + -характеристики при всех I > 0.

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

Описание формирования спирального течения со стоком от момента £ = 0 до достаточно больших значений времени возможно численно с помощью известного метода характеристик [8, 9]. Тем самым определится процесс поведения нестационарного течения, в том числе время его выхода на стационарный режим.

Пусть при некоторых значениях г0, гп известно решение стационарной задачи с°(г), п°(г), V°(г), где 0 < г0 < г < г^п.

Решение нестационарной задачи строится в областях С0, С1, представленных на рис. 4. При 0 < £ < ¿1 в области С0 решается задача (1), (3) с заданной функцией п0(£). Здесь ¿1 есть тот момент времени, в который рассчитанная радиальная скорость в точке г = гп становится равной (естественно, приближенно) значению п°(гп):

п(£,г ,т=Пп = п° (гш )•

Вычисленные к этому моменту времени £ = ¿1 при г0 < г < гп значения газодинамических параметров обозначаются как

и(1,г)\1=11 = и00 (г). (11)

При этом в расчетах при малых по модулю значениях щп величины с00 (гп), v00 (гп) получались достаточно близкими к единице и нулю соответственно.

После момента времени £ = ¿1 решение строится в полосе С1 (см. рис. 4):

: [г > £1; г0 < г < гп}.

Рис. 4. Область построения решения нестационарной задачи

На левой границе этой полосы С1, т. е. при т = т0, ставится условие, передающее непрерывный сток:

{п0(Ь), если п0(Ь) > п°(т0),

(12)

п°(т0), если п0(Ь) < п°(т0).

Иначе говоря, при достижении в точке т = т0 значения п°(т0) скорость стока полагается постоянной и равной скорости стока в стационарном течении.

На нижней границе полосы С\, т. е. при Ь = Ь1, задаются условия (11). На правой границе С\, при т = тп, задаются стационарные значения газодинамических параметров, определяемые значениями (11):

и(г,т%=Пп = и00 (тп). (13)

Система (1) имеет две звуковые характеристики С± и одну контактную С0:

С ± :

¿г 0 ¿г

— = и ± с, С : — = и, аЬ аЬ

вдоль каждой из которых вводится свое дифференцирование

(Й

с ±

д . . д а

с0

д_ д_

дЬ дт

и система (1) записывается в эквивалентном виде — в виде следующей системы обыкновенных дифференциальных уравнений:

ш

(Й

¿ь

(И

¿у

(И

с+

с-

¡1(т,С,п,Ь )\с+

с0

¡2(т,С,п,Ь )\

¡3(т,п,У )\с о

с-

где инварианты Римана Я, Ь задаются стандартным образом

2 2

г) = и + --—с, г) = и — --—с,

(7 - 1)

(7 - 1)

а правые части системы (14) имеют вид

(V \ си

/1 (г, с, и, V) = V [ —I-а)--,

тт

/2(г, с, и, у) = V + а^ + ¡3(г,и,ь) = -и + а^ .

Сп т

(14)

(15)

Определение газодинамических параметров по значениям инвариантов Римана осуществляется с помощью соотношений

с=^-Л(Н-Ь), и = \(П + Ь).

(16)

Для численного построения решения поставленной задачи применяется метод характеристик, который хорошо зарекомендовал себя при решении многих задач газовой динамики [8, 9].

Традиционный способ построения характеристической сетки, в узлах которой методом характеристик определяются значения искомых функций, приводит к построению треугольных ячеек с неравномерно расположенными узлами [8, 9]. Кроме этого, в случае системы, состоящей из более чем двух уравнений, необходимо применять интерполяцию [8,9], поскольку точка пересечения характеристики, имеющей промежуточный наклон, с линией, соответствующей предыдущему временному слою, в общем случае не совпадает с узлами характеристической сетки, ранее найденными по значениям максимального и минимального наклона характеристик.

В настоящей работе система (14) решается на прямоугольной сетке и поэтому используется стандартное для разностных схем обозначение И вектора значений искомых функций в точке (£ = £п ,г = гг), где

£п = пт, гг = г0 + гН.

На рис. 5 приведены две ячейки расчетной сетки.

Шаг сетки по пространственной переменной постоянен: Аг = Н. Выбор величины шага по временной переменной А£ = т определяется из неравенства

Н

т

\и\ + с'

(17)

выполнение которого гарантирует, что все три характеристики, выходящие из точки (£п+1, гг), пересекают прямую £ = £п в пределах отрезка [гг_1,гг+1 ] (см. рис. 5).

По значениям газодинамических параметров, заданных в точках (£п, гг_1), (£п, гг), (£п, гг+1), параметры газа в точке (£п+1 ,гг) определяются по следующему алгоритму. Вначале вычисляются величины А+, Л0 и А_ по формулам

А+ =

а(< + сп)

1 + а [(< + сп) - (<_ 1 + сп_ 1)]'

т

А0 =

—аиг

А_ =

1 + а(игп - ип_1)'

-а(ип + сп)

1 + а [(ип+1 - сп+1) - (игп - сгп)]

где а = —. С их помощью определяются значения Н

г+ = гг - А+ Н' г0 = гг + А°Н, ^

= гг + А Н

Рис. 5. Две ячейки расчетной сетки

координат точек, лежащих при £ = на отрезке [гг_ 1, гг+1] (см. рис. 5). Величины искомых функций в найденных точках с координатами (£?, г+), (£п, г0) и (£п, г_)

и(£п,г+) = Ип+, и(£п,г0) = И?о, и(£п,г_) = ип

вычисляются с помощью линейной интерполяции по значениям И™_^ И™ для точки (£?, г+) и Ип, И™+1 для точек (гп, г0), (гп, г_):

Ип+ = ип + А+ И?_! - Ип

И?0 = ип + А0 (И+1 - Ип

И?7 = ИП + А_ Иг+! - И").

Таким образом, три прямые, выходящие из точек (£?, г+), (£?, г0), (£?, ), с наклонами С+ -, С0- и С_-характеристик, вычисленными по значениям И"+, Ига0, ИТ соответственно,

Гг г Гг

проходят через точку (£П+1 ,гг) (см. рис. 5).

Далее с помощью разностной аппроксимации системы (14) находятся значения инвариантов Римана и функции v на следующем временном слое в точке (£п+1, гг):

К+1 = К+ + т/1 (г+, и'п+) ,

г ^ г /

Ь? +1 = + т/2 (г_ , и'п-) ,

г ^ г /

vгг+1 = Ф+т/з (г°, И^о),

г V г /

с помощью которых определяются

- ьг1), = £ (дг1 + ¿г1)

сГ1 = (дГ+1 _ ¿2+1), = I

Затем для улучшения точности расчета осуществляется пересчет значений искомых функций в точке (£П+1, гг) с переложением идеи традиционного для метода характеристик пересчета [8] на его данную модификацию. Для этого сначала вычисляются значения

д+= д[К+1 + сГ1) + (и? + С?)]

~0_ -а«+1+<)

2 + а(<+1 - и?)

д-= -д[(К+1 - сГ1) + (< - с?)} 2 + а[«+1-с-+1)-«-с-)]

по которым определяются уточненные значения

- ... _ Г+ ъ - ... _1_ Г0;

г+ = гг - А+Н, г° = гг + А°Н, г_ = г» + А Н,

и в этих точках (£?, Г+), (£?, Г0) и (£?, Гг ) вычисляются уточненные значения газодинамических параметров ( )

И+ = И? + А+ (и?_ 1 - Иг

и П = ип + л0 (Ип+1 - Ип

И- = ип + л - и+1 - ип

Далее определяются пересчитанные значения инвариантов Римана и функции V в точке (¿п+1, гг) по следующим формулам:

Затем находятся

К,

п+1

тп+1

г

¿г1 = Ц- + \ [л (гг, щ)+/2 {п, иг1);

«Г1 = е?» + ^ г 2

/э (К0, И "о) + /а (гг, и«+1);

с"+1 =

(7 - 1)

(нг1 - ¿Г1), «Г1 = ^ (К+1 + ¿Г1)

и производится переобозначение ИП+1 ^ ИП+1. Тем самым расчет значений газодинамических параметров в точке (£п+1, гг) заканчивается.

Расчет параметров газа на левой границе расчетной области при г = г0 (рис. 6) производится стандартным для метода характеристик способом. Значение и^+1 задано из условий (12). Вдоль С0-характеристики из точки (£ = 1п,г = г0) переносится значение функции V(Ь,г) и определяется значение v(n+1, вдоль С--характеристики из точки (£ = ,г = г-) переносится значение инварианта Ь(Ь,г) и, значит, становится известным значение ЬЩ+1. С использованием формулы (16) по заданным и^+1 и ЬЩ+1 сначала определяется значение Я^+1, а затем с£+1.

Для окончательного определения ИП также осуществляется соответствующий пересчет.

Для проверки правильности работы представленной программы, реализующей описанный алгоритм, численно строилось решение системы (14) с разными шагами по пространственной переменной К = 10-э, К = 10-4 при выполнении на каждом временном шаге условия (17). Расчеты показали определенную устойчивость и сходимость получаемых предложенным способом численных решений.

Если в качестве функций И00(г) (см. условия (11)) брались функции с =1, и = и0(г), V = 0, то в процессе счета происходил выход на стационарный режим и точное стационарное решение И0(г) в нестационарном расчете восстанавливалось с относительной погрешностью не более 1 %.

Рис. 6. Левая граничная расчетная ячейка

г 10 км 1 км 200 м 20 м(сток)

|"u|, км/ч 0.3 3.4 16.9 176.9

V, км/ч 0 6.7 33.9 353.7

Входные данные одного расчета нестационарного плоского спирального течения со стоком и выходом решения на стационарный режим в случае ф = п/6 следующие: h = 10-4, rin = 10 км, r0 = 20 м, \uo(rin)\ = 0.3 км/ч.

В таблице приведены размерные значения модуля радиальной и значения окружной составляющих вектора скорости газа в момент времени t = t2 = 6 ч 40 мин при различном размерном расстоянии. В рассматриваемом случае выход на стационарное течение наблюдается при t ~ 14 ч 10 мин, при этом значение окружной скорости в точке стока v(rin) = 654.5 км/ч.

Благодарим коллегу А.В. Боровских за полезные замечания, сделанные при выполнении данного исследования.

Список литературы

[1] Кочин Н.Е., Кибель И.А., Розе Н.В. Теоретическая гидромеханика. Ч. 1.М.: Физматгиз, 1963.

[2] Баутин С.П. Характеристическая задача Коши для квазилинейной аналитической системы// Диф. уравнения. 1976. Т. 12, № 11. С. 2052-2063.

[3] Баутин С.П. Математическая теория безударного сильного сжатия идеального газа. Новосибирск: Наука, 1997.

[4] Баутин С.П. Математическое моделирование сильного сжатия газа. Новосибирск: Наука, 2007.

[5] Баутин С.П. Характеристическая задача Коши и ее приложения в газовой динамике. Новосибирск: Наука, 2009.

[6] Петровский И.Г. Лекции об уравнениях с частными производными. М.: Физматгиз, 1961.

[7] Курант Р. Уравнения с частными производными. М.: Мир, 1964.

[8] Жуков А.И. Применение метода характеристик к численному решению одномерных задач газовой динамики // Тр. Математического ин-та им. В.А. Стеклова. М.: Изд-во АН СССР, 1960. Т. 58.

[9] Рождественский Б.Л., Яненко Н.Н. Системы квазилинейных уравнений и их приложения к газовой динамике. М.: Наука, 1968.

Поступила в редакцию 25 декабря 2009 г., с доработки — 23 августа 2010 г.