Научная статья на тему 'К вопросу о несимметрии и неединственности решения задачи об отрывном обтекании конической компоновки крыло корпус малого удлинения'

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

CC BY
198
44
i Надоели баннеры? Вы всегда можете отключить рекламу.
Журнал
Ученые записки ЦАГИ
ВАК
Область наук
Ключевые слова
НЕСИММЕТРИЯ / НЕЕДИНСТВЕННОСТЬ / ОТРЫВНОЕ ОБТЕКАНИЕ / КОНУС ТРЕУГОЛЬНОЕ КРЫЛО / МАЛОЕ УДЛИНЕНИЕ / РАСЧЕТЫ С УЧЕТОМ ВЯЗКОСТИ

Аннотация научной статьи по математике, автор научной работы — Воеводин Александр Владимирович

Проведены расчетные исследования отрывного обтекания компоновки круговой конус треугольное крыло малого удлинения в широком диапазоне определяющих параметров при симметричном расположении компоновки относительно набегающего потока. Использованы два подхода к решению задачи. В рамках теории тел малого удлинения получены неединственные симметричные и несимметричные решения. Уточнена область существования устойчивых решений разного типа. Проведены 3D-RANS CFD расчеты компоновки малого удлинения с носовой частью в виде кругового конуса с треугольным крылом с учетом вязкости. Показано, что в этом случае задача имеет единственное решение, которое симметрично при малых углах атаки и несимметрично при больших. При относительных углах атаки меньше 2 коэффициенты подъемной силы, полученные с помощью двух подходов, очень близки, что дает оценку области применимости теории тонких тел. Приведены картины течения в сечении конической части компоновки, иллюстрирующие детали обтекания.

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Текст научной работы на тему «К вопросу о несимметрии и неединственности решения задачи об отрывном обтекании конической компоновки крыло корпус малого удлинения»

Том X Ь

УЧЕНЫЕ ЗАПИСКИ ЦАГИ 2009

№ 6

УДК 629.735.33.015.3:533.695.12

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

А. В. ВОЕВОДИН

Проведены расчетные исследования отрывного обтекания компоновки круговой конус — треугольное крыло малого удлинения в широком диапазоне определяющих параметров при симметричном расположении компоновки относительно набегающего потока. Использованы два подхода к решению задачи. В рамках теории тел малого удлинения получены неединственные симметричные и несимметричные решения. Уточнена область существования устойчивых решений разного типа. Проведены 3D-RANS CFD расчеты компоновки малого удлинения с носовой частью в виде кругового конуса с треугольным крылом с учетом вязкости. Показано, что в этом случае задача имеет единственное решение, которое симметрично при малых углах атаки и несимметрично при больших. При относительных углах атаки меньше 2 коэффициенты подъемной силы, полученные с помощью двух подходов, очень близки, что дает оценку области применимости теории тонких тел. Приведены картины течения в сечении конической части компоновки, иллюстрирующие детали обтекания.

Ключевые слова: несимметрия, неединственность, отрывное обтекание, конус — треугольное крыло, малое удлинение, расчеты с учетом вязкости.

Исследованию невязкого стационарного отрывного обтекания конической компоновки крыло — корпус малого удлинения посвящено много работ. При условии малости угла атаки и угла при вершине треугольного крыла исходная трехмерная стационарная задача обычно сводится к двумерной нестационарной, и затем ищется ее автомодельное решение. В первых работах предлагались и усовершенствовались различные модели вихревой пелены, сходящей с кромки крыла [1—4]. В наиболее полной из них внешняя часть пелены моделировалась конечным числом отрезков с распределенной на них завихренностью, а внутренняя часть вихревой спирали заменялась дискретным вихрем. Для определения формы и интенсивности пелены использовалось уравнение ее эволюции [5]. Существенным моментом большинства работ было использование конформного отображения контура тела в поперечном сечении на отрезок, расположенный вдоль набегающего потока, или на окружность. Это позволяло автоматически выполнить условие непротекания. Можно также отметить работу [6], где на контуре тела распределялись аэродинамические особенности, интенсивность которых определялась из условия непротекания. Такой подход позволяет исследовать конические компоновки, контур поперечного сечения которых не может быть отображен на отрезок или окружность с помощью простых формул. Для решения полученных уравнений использовался итерационный метод или метод установления [7—9].

С помощью описанной методики были проведены многочисленные расчеты, например, компоновки треугольного крыла с коническим фюзеляжем, имеющим в сечении полукруг [10], а также крыла, имеющего в сечении дужку окружности [11]. В работе [12] приведены результаты расчетов обтекания треугольного крыла при наличии скольжения.

При исследовании симметричного обтекания компоновки круговой конус — треугольное крыло [13] была обнаружена неединственность решения задачи в определенном диапазоне угла

атаки и отношения диаметра фюзеляжа к размаху крыла. В работе [14] было показано, что если не требовать заранее симметричности решения (т. е. решать задачу для полной компоновки, а не для ее половины), то при симметричном расположении компоновки относительно набегающего потока задача имеет как симметричные, так и несимметричные решения, среди которых имеются неустойчивые. Однако исследование устойчивости проведено для простейшей модели вихревой пелены («вихрь — разрез»), поскольку в этом случае такое исследование можно выполнить в аналитическом виде.

Первая часть настоящей работы посвящена исследованию в рамках теории тонких тел области существования устойчивых симметричных и несимметричных решений задачи об отрывном обтекании компоновки круговой конус — треугольное крыло при симметричном расположении компоновки относительно набегающего потока. Во второй части проведены расчеты компоновки конечной длины с носовой частью в виде конуса с треугольным крылом методом 3D-RANS.

Постановка задачи. Исследуется отрывное обтекание компоновки круговой конус — треугольное крыло малого удлинения. Жидкость будем считать невязкой и несжимаемой, сдвиговые слои — вырожденными в поверхности тангенциального разрыва скорости, а линии отрыва — прикрепленными к острой кромке крыла. Полуугол при вершине крыла 5 и угол атаки а малы: 5 ~ а = о(1). Ось X системы координат OXYZ направлена по оси компоновки, ось Y — вверх, и

ось Z лежит в плоскости крыла. Согласно стандартной процедуре исходная трехмерная задача (в переменных X, Y, Z) сводится к двумерной нестационарной об отрывном обтекании равномерно расширяющегося тела, представляющего собой поперечное сечение компоновки (в переменных t = X/U„ , Y, Z, где U„ — скорость набегающего потока). Далее вводятся автомодельные переменные х и у, так что Y = yU^ttg 5 и Z = xUjttg 5, и комплексная переменная z = х + iy. В плоскости z скорость набегающего потока (относительный угол атаки) а0 = sin а/tg 5 = а/5, кромки крыла находятся в точках —1 и 1, а радиус фюзеляжа обозначим через m <1 (рис. 1).

Тогда, согласно [4, 5], уравнение эволюции вихревой пелены, которое является следствием условий отсутствия скачка нормальной к пелене составляющей скорости и отсутствия скачка давления поперек пелены, записывается в виде:

Рис. 1. Контур тела и вихревая пелена в плоскости z

_ . аг dw

г -Аф—-------- =—,

d (Аф) dz

где w(г) — безразмерный комплексный потенциал течения, ф = Яе(w), а Аф = Г есть скачок вещественной части комплексного потенциала поперек пелены. В точке (-1, 0) Аф = Gl, а в точке (1, 0) Аф = в2, где Gl и в2 — полные циркуляции левой и правой вихревых пелен.

Для построения комплексного потенциала течения отобразим внешность контура тела на внешность окружности единичного радиуса в плоскости комплексной переменной ц = £, + /'П. Такое преобразование задается формулой:

1 + т

2 (

г =-

+

1 + т2 ( 1 і

Ц + -

т

Учитывая, что на бесконечности ^^ =—и используя метод отражений, комплексный по-

dz 1 + т2

тенциал запишем в виде:

w

(ц) = т2 1п г -/а(

т2 +1 (

1

Ц —

ц

л

+

— Г 1п Ц Ц(1Cp)dАф + — Г 1пГ Т/ dАф. (2)

2п і* 1 2п/і 1

ц-ц(Аф)

С! ц-

ц(Аф)

2П ‘Сг ц- 1

ц(Аф)

Построенный таким образом потенциал удовлетворяет условию непротекания на теле и на бесконечности ведет себя как -/ао г.

Учитывая (2) и вводя новую переменную х = 1 ^■Аф, уравнение (1) приводится к виду:

G

2 а ( г і т2 . т2 +1 ^

(1 -Х)2^-1^-1 = —— /а0----

'• ’ ах 1 1 -X У г 0 2

1+4

ц ,

Л

а ц

+

+

2п /

1 ( 1 1 і 1 ( 1 1 і

ц і хх 1 а х+G7 Г ц-ц2 (х) ,, 1 а х

12 0 ц(х) у і 0 1 ц2(Х)у

ац

(3)

Индексы «1» и «2» в подынтегральных выражениях указывают, что интегрирование ведется по контуру левой или правой пелены.

Условие Чаплыгина — Жуковского на острых кромках крыла, где производная йц/йх обращается в бесконечность, дает еще два вещественных уравнения, из которых при известной геометрии пелен можно получить интенсивности 01 и G2:

I 2 \ 1 (1 -<а (^)-П12 (X)) 1 (1 -(£ (х)-п2 (X))

2па0 (т +1) + 01 Г------------------ йХ + 02 Г---— йХ = 0, (4)

о (-& (х)) +П2 (х) о (-^2 (х)) + П (х)

2 ( 2 1) 0 1 (1 -^2(Х)-П2(х)) 0 1 I1 -^2(Х)-П2(х)) о (5)

2па0 (т +11 + 011----------■------- йХ + 02 I------------------------------------■-— йХ = 0. (5)

о ( + ^1 (х)) +Л2 (х) 0 (1 + ^2 (х)) +п2 (х)

Метод численного решения. Для численного решения задачи каждая пелена разделялась на внешнюю часть, которая разбивалась на конечное число отрезков Л и N2, и внутреннюю, которая заменялась дискретным вихрем (ядро). Конец последнего отрезка внешней части соеди-

нялся разрезом с дискретным вихрем для преодоления неоднозначности потенциала скорости (показаны пунктиром на рис. 1). В соответствии с этим интегрирование в уравнениях (3) — (5) проводилось соответственно от нуля до хи от нуля до хд, а интегралы по внутренним частям пелен заменялись соответствующим воздействием от дискретных вихрей. Интегралы по отрезкам пелены, непосредственно примыкающим к кромкам (1), вычислялись с использованием асимптотики вблизи острой кромки тонкого крыла, согласно которой в рассматриваемой задаче

ц (х) = — 1 — й^х1^2 + //^х +...,

М"2 (х)=1+й^х1^2+\Ь2х+...

Вещественные константы а-^, а2, Ь1 и Ь2 определялись численно.

Для определения координат ядер использовалось стандартное условие равенства нулю полной силы, действующей на систему вихрь — разрез, которое имеет следующий вид:

2хС — = й- , (6)

йх с

где хс есть координаты ядра 1 или 2, а хд — координаты конца внешней части соответствующей пелены. Комплексно сопряженная скорость вычислялась в соответствующей точке Хс .

Уравнения (3) — (6) решались итерационным методом. При этом уравнение (3) выполнялось в средних точках отрезков пелен, а интегрирование проводилось методом трапеций, имеющим второй порядок аппроксимации. В качестве начального приближения использовалась модель пелены, в которой внешняя часть состояла только из асимптотического, примыкающего к кромке участка. Решение для такой модели быстро сходится, и далее количество отрезков внешней части и ее протяженность последовательно наращивались до формирования примерно одного витка спирали. После этого проводилось вычисление сил, действующих на компоновку. Более подробно детали конечно-разностного представления уравнений и алгоритма проведения расчетов приведены в работе [13].

Результаты расчетов. Как уже упоминалось ранее, при симметричной постановке задачи в рамках теории тонких тел решение получается неединственным. В определенном диапазоне углов атаки и относительных радиусов фюзеляжа имеются два устойчивых и одно неустойчивое к симметричным возмущениям решения. При исследовании решений в несимметричной постановке [14] показано, что одно из симметричных решений неустойчиво к несимметричным возмущениям. Кроме того, имеется набор устойчивых и неустойчивых несимметричных решений. Поскольку эти исследования проводились в рамках упрощенной модели вихревой пелены (вихрь — разрез), остается невыясненным вопрос — как изменится ситуация при применении более полной модели. В настоящей работе проведены расчеты течения в широком диапазоне параметров задачи для полной модели пелены (пелена — вихрь — разрез) в несимметричной постановке.

Рассмотрим результаты этих исследований. На рис. 2 и 3 показаны зависимости коэффициентов подъемной и боковой силы от относительного угла атаки а0 для разных значений радиуса

фюзеляжа т ( = tg 5). Индекс «5» соответствует симметричному решению, а «ш» — несимметричному. Устойчивое симметричное решение существует для всех значений т и а0 <а01 (т).

При больших а0 симметричное решение теряет устойчивость. В этом случае при т < т* ~ 0.8 решение становится неустойчивым к несимметричным возмущениям, но в их отсутствии может быть продолжено на большие углы атаки. Если т > т* и а0 >а01 (т), симметричное решение

неустойчиво как к несимметричным, так и к симметричным возмущениям.

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

Рис. 2. Зависимости коэффициента подъемной силы от относительного угла

атаки

200

150

100

50

С=/К2

-о- т =0.3 -о- т — 0.4 т = 0.5 — т = 0.7 т =0.75 -с- т =0.9

^5

0 2 4 6 ^0

Рис. 3. Зависимости коэффициента боковой силы от относительного угла атаки

Кривая а01 (т) на этой диаграмме качественно повторяет аналогичную зависимость, приведенную в работе [14], однако, как и предполагалось, она оказалась смещенной вниз по углам атаки. Это отличие вызвано тем, что в [14] использовалась упрощенная модель вихревой пелены. Кроме того, на рис. 4 нанесена кривая а02 (т), соответствующая появлению несимметрии, и

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

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

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

а0 ч

\ '© у г—ЛЮ

©

/ а01(?и) ©

О а 02 («0

О 0.2 0.4 0.6 0.8 1

Рис. 4. Области существования устойчивых решений задачи:

□ — появление несимметричного решения;-------потеря устойчивости сим-

метричного решения; о — несимметричное решение 1; А — несимметричное

решение 2

Рис. 5. Конфигурация вихревой пелены:

а — т = 0.8, а0 = 3.5; б - т = 0.8, а0 = 5.4; в - т = 0.98, а0 = 5.2

-------симметричное решение

А — решение 1

□ — решение 2

На рис. 5 показана форма вихревых пелен для некоторых значений параметров а0 и т, соответствующих областям 3, 4 и 5. Видно, что конфигурации пелен, соответствующих разным решениям, существенно отличаются.

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

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

Обе компоновки имеют одинаковый фюзеляж: длина конической части 1 м, длина цилиндрической части 0.5 м, а хвостовая часть представляет собой полусферу. Крыло первой из них (рис. 6) имеет полу-угол при вершине крыла 8 = 5°, при этом отношение диаметра фюзеляжа к размаху крыла т=0.9. Полу-угол при вершине второго крыла 8 = 6°, что соответствует параметру т = 0.75. Так как существенного качественного отличия в результатах расчетов для Рис. 6. Геометрия компоновки этих двух вариантов не обнаружено, далее в основ-

ном приводятся данные для т= 0.9. Расчетная область имела размеры: по X от -1.5 до 5 м, по У от -1.3 до 1.7 м и по 2 от -1.5 до 1.5 м и содержала около 1 млн. ячеек. Носик фюзеляжа располагался в начале координат.

Решались осредненные по Рейнольдсу уравнения Навье — Стокса с ББТ моделью турбулентности. Набегающий поток имел скорость 50 м/с, плотность 1.2 кг/м3 и вязкость

1.83 10-5 кг/м • с. Расчеты проводились для углов атаки до 45°, что соответствует а0 = 8.082. На рис. 7 и 8 показаны полученные зависимости коэффициентов подъемной и боковой сил в сечении Х= 0.7 м от относительного угла атаки при т = 0.9. Для удобства сравнения с результатами теории тонких тел (бЫ) их величины поделены на К . На рис. 7 пунктиром нанесена кривая, соответствующая расчету половины компоновки (симметричная постановка задачи). Видно, что коэффициенты подъемной силы в двух рассмотренных случаях очень близки. Отличие наблюдается только в диапазоне а = 26 -^37°, который соответствует максимальной по модулю боковой силе в несимметричном решении.

Во всем исследованном диапазоне угла атаки реализуется единственное решение. При этом до угла атаки 12° (0 = 2.38) оно симметрично и коэффициенты подъемной силы, полученные

по теории тонких тел и методом CFD, практически совпадают. На больших углах атаки решение теряет симметрию, и появляется боковая сила. Характер ее поведения соответствует экспериментальным данным [15, 16]. Какие-либо перемежающиеся вихревые структуры, которые наблюдаются при обтекании тонких тел вращения, не обнаружены. На рис. 9 показаны линии тока, сходящие с кромок крыла на конической части компоновки при а = 30°, m = 0.9. Здесь хорошо видно, что течение несимметрично и в этой области даже при таком большом угле атаки близко к коническому. Известно, что в конических течениях линии тока в окрестности ядра вихря практически цилиндрические. На рис. 9 можно заметить небольшое отклонение этих линий тока от цилиндрической поверхности (уширение), что может быть вызвано большим углом атаки, способом интегрирования линий тока и их большой протяженностью.

На рис. 10 представлены картины полей завихренности для компоновки с m = 0.75 в поперечном сечении Х= 0.7 м для углов атаки 10 и 15°.

Здесь точками изображены вихревые пелены, рассчитанные по теории тонких тел. При малых углах атаки оба использованных метода дают очень близкие структуры течения (как и интегральные характеристики). При а = 15° в CFD расчете основная область завихренности уходит выше по потоку, и становятся заметными несимметрия и вторичный отрыв с поверхности конуса. При еще больших углах атаки отличия картин течения усиливаются.

Для исследования картин течения в сечении X = const сформируем следующее модифицированное поле скоростей:

Рис. 9. Линии тока, сходящие с кромок треугольного крыла (а = 30°, т = 0.9)

U1 = U -

V = V-

-U„ cos а R

Y

Y zz2)x’

(y 2 +z2 X’

где и, V и Ж — компоненты скорости; Я — радиус фюзеляжа в сечении. На рис. 11 показаны «линии тока», соответствующие такому полю скоростей для т =0.9, Х= 0.7 м и разных углов атаки.

Рис. 10. Поля завихренности в плоскости X = 0.7 м (m= 0.75)

Рис. 11. Картины течения в плоскости Х = 0.7 м (т = 0.9)

В диапазоне а = 9 + 12° на верхней поверхности фюзеляжа имеется вторая заметная область отрыва. Направления движения газа в этой области и в области основного отрыва совпадают. С дальнейшим увеличением угла атаки течение становится несимметричным, второй отрыв исчезает, и появляется вторичный отрыв с фюзеляжа, расположенный ближе к крылу. Кривая су (а)

в окрестности а = 12° существенно меняет наклон. При а « 25° течение симметризуется и боковая сила пропадает. Далее боковая сила меняет знак и при а « 32° достигает минимума. Еще один ее переход через ноль происходит при а~ 38.5°.

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

Проведены расчеты компоновки конечной длины методом 3Б ЯЛК8. Показано, что решение при больших углах атаки с учетом вязкости имеет существенно иной характер. Решение получается единственным, при этом на малых углах атаки оно симметрично, коэффициент подъем-

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

ЛИТЕРАТУРА

1. Legendre R. Ecoulement au voisinage de la point avant d’une aile a forte fleche aux incidences moynnes // Rech. Aeron. 1952. N 30.

2. Brown C. E., Michael W. H. Effect of leading-edge separation on the lift of a delta wing // JAS. 1954. V. 21(2).

3. Жигулев В. Н. Некоторые задачи нелинейной теории крыла // Труды ЦАГИ.

1955, Ж-688.

4. Smith J. H. B. Improved calculation of leading-edge separation from slender delta wings // Royal Aircraft Establishment. Tech. Rep., 1966, 66070.

5. НикольскийА. А. О «второй» форме движения идеальной жидкости около обтекаемого тела (исследование отрывных вихревых потоков) // Доклады АН СССР. 1957.

Т. 116, № 2.

6. Воеводин А. В., Судаков Г. Г. Проекционный метод расчета характеристик отрывного обтекания тел идеальной жидкостью // Ученые записки ЦАГИ. 1986. Т. XVII, № 5.

7. Sacks A. H., Lundberg R. E., Hanson C. W. A theoretical investigation of the aerodynamics of slender wing-body combination exhibiting leading-edge separation // NASA.

1967, CR-719.

8. Молчанов В. Ф. О реализации метода плоских сечений в нелинейной теории крыла // Ученые записки ЦАГИ. 1974. Т. V, № 2.

9. СудаковГ. Г. Расчет отрывного течения около тонкого треугольного крыла малого удлинения // Ученые записки ЦАГИ. 1974. Т. V, № 2.

10. ВоеводинА. В. Отрывное обтекание конической комбинации крыла малого удлинения с несимметричным фюзеляжем // Ученые записки ЦАГИ. 1983. Т. XIV, № 4.

11. В оеводинА. В. О влиянии формы в плане и кривизны поверхности крыла малого удлинения на характеристики его отрывного обтекания // Труды ЦАГИ. 1986, вып. 2317.

12. P u 11 i n D. I. Calculations of steady conical flow past a yawed slender delta wing with leading-edge separation // ARC R & M. 1975. N 3767.

13. Воеводин А. В. Исследование неединственности решения задачи об отрывном обтекании системы крыло-фюзеляж малого удлинения // Ученые записки ЦАГИ. 1979.

Т. X, № 1.

14. ГоманМ. Г., Захаров С. Б., Храбров А. Н. Симметричное и несимметричное отрывное обтекание крыла малого удлинения с фюзеляжем // Ученые записки ЦАГИ.

1985. Т. XVI, № 6.

15. Andrew B., Wardlaw Jr. Multivortex model of asymmetrical shedding on slender bodies at high angle of attack // AIAA Paper. 1975, 75-123.

16. Захаров С. Б., Зубцов А. В. Экспериментальное исследование отрывного обтекания треугольного крыла малого удлинения // Ученые записки ЦАГИ. 1988. Т. XIX, № 1.

Рукопись поступила 11/II2009 г.

i Надоели баннеры? Вы всегда можете отключить рекламу.