Научная статья на тему 'Влияние установившегося вращения на аэродинамические характеристики профиля при наличии отрыва потока'

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

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

Аннотация научной статьи по механике и машиностроению, автор научной работы — Волков А. В., Ляпунов С. В., Храбров А. Н.

С использованием приближенного аналитического и численных подходов для широкого диапазона углов атаки исследуются нелинейные аэродинамические характеристики профиля NACA 0012 при числе Маха М=0,15 и числе Рейнольдса Re=1·106 при его установившемся вращении по тангажу с различными угловыми скоростями относительно произвольной точки хорды. Показано, что установившееся вращение оказывает заметное влияние на коэффициент подъемной силы и момента тангажа профиля. Вращение профиля по тангажу с положительной угловой скоростью относительно четверти хорды приводит к более раннему развитию отрыва потока по сравнению с профилем без вращения. Вращение около трех четвертей хорды приводит к обратному эффекту. Результаты расчетов с использованием различных методик хорошо согласуются друг с другом.

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

Похожие темы научных работ по механике и машиностроению , автор научной работы — Волков А. В., Ляпунов С. В., Храбров А. Н.

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

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

Том XXXIV

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

№ 3—4

УДК 629.735.33.015.3.025.73:532.526

ВЛИЯНИЕ УСТАНОВИВШЕГОСЯ ВРАЩЕНИЯ НА АЭРОДИНАМИЧЕСКИЕ ХАРАКТЕРИСТИКИ ПРОФИЛЯ ПРИ НАЛИЧИИ ОТРЫВА ПОТОКА

А. В. ВОЛКОВ, С. В. ЛЯПУНОВ, А. Н. ХРАБРОВ

С использованием приближенного аналитического и численных подходов для широкого диапазона углов атаки исследуются нелинейные аэродинамические характеристики профиля КАСЛ 0012 при числе Маха М=0,15 и числе Рейнольдса Яе = 1 • 106 при его установившемся вращении по тангажу с различными угловыми скоростями относительно произвольной точки хорды. Показано, что установившееся вращение оказывает заметное влияние на коэффициент подъемной силы и момента тангажа профиля. Вращение профиля по тангажу с положительной угловой скоростью относительно четверти хорды приводит к более раннему развитию отрыва потока по сравнению с профилем без вращения. Вращение около трех четвертей хорды приводит к обратному эффекту. Результаты расчетов с использованием различных методик хорошо согласуются друг с другом.

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

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

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

поверхности совпадает с задней кромкой профиля) переходят в известные формулы теории тонкого профиля [3].

Для численного исследования отрывного обтекания профиля при наличии установившегося вращения по тангажу использовались два различных подхода с соответствующими программами расчета. Одна программа была разработана в ЦАГИ и широко используется для расчета отрывного обтекания профиля, в том числе со щелевой взлетно-посадочной механизацией [4]. Для сравнения была выбрана программа Хбэй [5], [6]. Эта программа разработана в Массачусетсском технологическом институте (США) для расчета отрывного обтекания односвязных профилей при дозвуковых скоростях полета. Она получила широкое распространение за рубежом благодаря тому, что распространяется в сети Интернет в качестве бесплатного математического обеспечения.

При адаптации этих программ для расчета обтекания вращающихся профилей был применен так называемый метод искривленных моделей. Идея этого метода принадлежит профессору В. П. Ветчинкину [7]. Она заключается в том, что для определения аэродинамических нагрузок на вращающемся профиле при испытаниях в прямолинейном потоке аэродинамической трубы используется специальным образом изготовленная искривленная модель. Искривление модели производится так, чтобы местные углы атаки в соответствующих точках неискривленной модели, движущейся с заданной угловой скоростью, и искривленной модели в прямолинейном потоке совпадали. Строгое обоснование данного подхода в линейном случае рассматривается в работе [8]. В нелинейном случае этот метод является приближенным. Экспериментальное подтверждение данного метода рассматривается в работе [9], где проведено сравнение данных, полученных в широком диапазоне углов атаки а и угловых скоростей тангажа О для модели прямоугольного крыла малого удлинения X = 0,5 с профилем КЛСЛ 0012 по методу искривленнных моделей в аэродинамической трубе ЦАГИ, с данными эксперимента для неискривленной модели той же формы, полученными на ротативной машине ЦНИИ им. А. Н. Крылова. Оба подхода дают весьма близкие результаты как для безотрывных, так и для отрывных режимов обтекания.

С использованием этих подходов для широкого диапазона углов атаки в статье исследуются нелинейные аэродинамические характеристики профиля КЛСЛ 0012 при числе М = 0,15 и числе Яе = 1- 106 при его установившемся вращении по тангажу с различными угловыми скоростями относительно точки, имеющей различную координату вдоль хорды.

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

скорость изменения угла атаки. В дальнейшем было бы интересно сравнить полученные в статье приближенные результаты с решениями уравнений Навье — Стокса.

1. Рассмотрим кратко различные подходы для расчета нелинейных аэродинамических характеристик профиля при наличии установившегося вращения по тангажу. Остановимся сначала на использовании теории тонкого профиля (линейной теории кавитации) с учетом отрыва потока в виде зоны постоянного давления Кирхгофа. Схема течения показана в верхней части рис. 1. Профиль с хордой 6 = 1

■го

Ср/а

10

20

0 к і і і і і і

\ *1=0,3

1 Г Г 1 | Г Г Г 1 :

Рис. 1. Учет отрыва потока в теории тонкого профиля

обтекается потоком со скоростью Уо = 1, направленной под углом атаки а. Профиль

может вращаться по тангажу с постоянной безразмерной угловой скоростью О вокруг точки хорды с продольной относительной координатой х0. При этом, конечно, траектория движения профиля не является прямолинейной. На верхней поверхности профиля в точке с относительной координатой х5 начинается отрыв потока, моделируемый полубесконечной зоной постоянного давления. В качестве упрощающих предположений принимается, что зона отрыва узкая (вследствие чего граничные условия могут быть снесены на ось, а не удовлетворяться на поверхности неизвестной формы) и след за профилем прямолинейный. При этих предположениях может быть получена [1], [2] аналитическая формула для распределения давления по поверхности тонкого симметричного профиля в зависимости от нормальной скорости на части его поверхности, обтекаемой безотрывно,

и координаты начала отрываа потока на верхней поверхности хя. В нижней части рис. 1, например, показано распределение давления по профилю, получаемое в такой задаче для заданных значений х5 = 0,3 и О = 0. Видна особенность на передней кромке профиля, присущая

теории тонкого профиля, и постоянное распределение давления на верхней поверхности за точкой отрыва (коэффициент давления равен нулю вследствие полубесконечности зоны Кирхгофа). Для пластинки это распределение давления может быть проинтегрировано аналитически, что приводит к явным формулам для аэродинамических коэффициентов подъемной силы су и момента тангажа т2 относительно условного центра тяжести Хо. Формула для су

имеет следующий вид:

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

Если в этих формулах положить х5 = 1, что соответствует безотрывному обтеканию профиля (начало отрыва потока на задней кромке), можно получить следующие выражения:

Эти выражения известны из линейной безотрывной теории тонкого профиля [3]. Из них, например, следует, что коэффициент подъемной силы при безотрывном обтекании не зависит от угловой скорости вращения профиля при Хо = 3/4 . Если же профиль вращается относительно

уп (х) = -а + О(х - х0 )

(1)

ш.

паа (1+^) (5 -6^+5х* - 16хо )-128 (1+^) (25 - 80х0+

+64х^ - 44^х^ + 96х0^/Х^ + 54х!, - 80хохх - 44х^^Х^ + 25х2)

(3)

середины хорды, то при малых углах атаки от угловой скорости вращения не зависит коэффициент

т.

Таким образом, с помощью формул (1) и (2) можно приближенно рассчитывать аэродинамические нагрузки на профиле при произвольных значениях параметров а, О и Хо, если задано значение положения точки отрыва потока на его верхней поверхности х*,, которое в свою очередь зависит от многих параметров, включая и три предыдущих. В дальнейшем в настоящей работе проводится сравнение результатов расчетов аэродинамических нагрузок по данным формулам с численными результатами, полученными по программам ЦАГИ и ХАзП. Зависимость х*, (а, О) при этом бралась из результатов численных расчетов.

По аналитическим формулам (1) и (2) можно вычислить также и соответствующие вращательные производные с°, т° . При О ^ 0 для производной коэффициента подъемной силы

имеем

—С- ^ + Л(1 + ПТ)2 (15 - 18хОЛ/х7 + 15х, - 16х,)

dО 2 —О 32' ^ ^ * * *>

и для коэффициента момента тангажа

—тг па (1 + >/ хя)/ I— \—** п / /—\2(~. „„ 2

(4)

- 80х0 + 64х0 -

—О 4 ^хх v 44 '—О 128^ * V и и (5)

-44д/х^ + 96х0л1х~ + 54х* -80х0х* -44х*^/х~ + 25х2).

Из формул (4) и (5) видно, что вращательные аэродинамические производные существенным образом зависят от функции .

2. Как уже упоминалось во введении, для расчета течения около профиля КЛСЛ 0012 при наличии установившегося вращения с учетом вязкого турбулентного пограничного слоя и возможного отрыва потока в работе используются две программы. Обе программы (ЦАГИ и ХйэП) для решения задачи используют зональный подход, который предполагает, что все вязкие силы действуют только в тонком пристеночном слое. Это позволяет разбить расчетную область на две части — невязкую, где решаются уравнения движения идеального газа, и вязкую, где уравнения движения вязкого газа значительно упрощаются благодаря тонкости соответствующего слоя. Согласование решений во внешней (невязкой) и внутренней (вязкой) областях достигается итерационно с учетом взаимодействия этих областей. Рассматриваемые программы позволяют в рамках зонального подхода проводить расчеты обтекания профилей в широком диапазоне углов атаки, чисел М и Яе, включая случаи развитых диффузорных отрывов.

В программе, разработанной в ЦАГИ, внешняя невязкая задача решается конечноразностным методом для уравнения потенциала. При этом используется конформное преобразование, переводящее внешность профиля во внутренность круга. Данный метод позволяет быстро получать решения, включая режимы с образованием местной сверхзвуковой зоны. Уравнения турбулентного пограничного слоя с учетом соответствующей полуэмпирической методики перехода от ламинарного к турбулентному течению решаются интегральным методом. Согласование решений для невязкого и вязкого течений производится на основе метода вязко-невязкого взаимодействия. Подробное описание методики расчета приведено в работе [4]. Результаты расчетов по данной методике обычно демонстрируют удовлетворительное согласование с экспериментальными данными в широком диапазоне изменяемых параметров.

В программе ХйэП внешнее невязкое течение находится на основе простого панельного метода. Вдоль верхней и нижней поверхностей профиля располагаются панели с линейным

изменением завихренности. Сжимаемость потока учитывается введением поправки Кармана — Цзяна, что не позволяет в этой программе рассчитывать течение с образованием местных сверхзвуковых зон. Пограничный слой также рассчитывается интегральным методом с

использованием вм критерия перехода. Этот критерий применим для предсказания турбулентного перехода, если доминирующим механизмом развития турбулентности является экспоненциальный рост двумерных волн Толлмина — Шлихтинга вследствие потери их устойчивости в линейном приближении. Этот эффект наблюдается для большинства профильных течений. Для применения данного метода пользователь должен задать параметр ^сг;1, который является логарифмом коэффициента усиления наиболее возрастающей частоты. Этот параметр зависит от степени турбулентности потока, набегающего на профиль. Для стандартной аэродинамической трубы авторы программы рекомендуют выбирать ^сг;1 = 9 (АДТ с высоким уровнем турбулентности — ^сг;1 = 4, натурное течение около крыла планера — ^с1.;1 = 12 -14).

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

Обе рассматриваемые методики численного расчета были разработаны без учета установившегося вращения профиля. Для включения этого эффекта рассмотрим граничные условия на поверхности профиля при его установившемся вращении по тангажу. В связанной системе координат нормальная компонента скорости на поверхности вращающегося тонкого профиля должна иметь следующий вид [3]:

^ =^х — -V,----------О

п х У

/ ч —х —У

(х - х0 )— + У— 0

(6)

где Ух = V) соб а и V, = V бш а — горизонтальная и вертикальная компоненты скорости движения профиля V), О — угловая скорость вращения по тангажу, у(*) и х(*) — параметрически заданная форма профиля.

Будем рассматривать только вращения с малой угловой скоростью О, тогда ее влияние может быть смоделировано малыми искривлениями поверхности исходного профиля. Будем искать уравнение поверхности искривленного профиля в виде

у' = У(*) + §(*) е х' = х(*) .

Так как искривленный профиль движется прямолинейно, нормальная скорость на его поверхности определяется выражением

^ = гЛ- Уу—=гх ( ^—*1-Гу^х. (7)

п х 1 У йс х йс йс У 7 К '

Потребуем выполнения условия равенства нормальных скоростей в соответственных точках. Тогда из сравнения выражений (6) и (7) будем иметь

— = -О

, чах ау

(х - х0 )— + У— V 0 / ^

(8)

что, очевидно, приводит к следующему виду функции, задающей искривление поверхности:

/ \2 2 ( х - х0 ) + У

8 = -О^-------01-------------------------------------------—. (9)

2

Выражение (9) совпадает в линейном приближении при у << 1 с выражением для искривленной поверхности профиля, полученным для линейного случая в работе [8]. При изготовлении экспериментальных моделей среднюю линию изгибают также по параболе.

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

поверхностью, деформированной в соответствии с выражением (9). При этом граничные условия на профиле будут удовлетворены с точностью О(О).

Неучтенным, в частности, остается эффект искривления следа за профилем.

3. Обратимся к обсуждению результатов расчетных исследований. В работе численно исследовалось обтекание вязким

турбулентным потоком профиля КЛСЛ 0012 при различных угловых скоростях его вращения по тангажу. При изучении влияния вращения важно учитывать не только величину угловой скорости, но и положение координаты центра вращения. Рассмотрены три случая положения мгновенного центра вращения профиля х0 = 0,25; 0,5 е 0,75 . Для каждого положения х0 рассчитаны вращения с безразмерной угловой скоростью тангажа О = 0; ± 0,05 и ± 0,1. Для каждого случая с учетом известного аналитического выражения для формы профиля КЛСЛ 0012 и выражения (9) была рассчитана деформированная форма изогнутого профиля. Таким образом, расчеты проводились для двенадцати изогнутых профилей и одного исходного симметричного. Все расчеты выполнены для чисел М = 0,15 и Яе = 1-106, соответствующих стандартным условиям эксперимента в АДТ малых дозвуковых скоростей.

На рис. 2 показано сравнение результатов, полученных различными

Рис. 2. Сравнение результатов расчета аэродинамических характеристик невращающегося профиля ЫЛСЛ 0012 при М=0,15 и Яе=106 при наличии отрыва потока с помощью различных подходов

методами для

и

исходного

невращающегося профиля. Условный центр тяжести в этих расчетах находился в точке

х0 = 0, 25 ,

относительно

этой

точки

Рис. 3. Влияние установившегося вращения около точки х0 = 0,25 на аэродинамические характеристики профиля по программам ЦАГИ и ХМ1

рассчитывался коэффициент момента тангажа. В расчетах по линеаризованной теории тонкого профиля с использованием формул (1) и (2) зависимость координаты начала отрыва потока на верхней поверхности от угла атаки хх (а) бралась по результатам расчета

пограничного слоя в методике ЦАГИ. Видно, что все три метода дают результаты, очень близкие между собой.

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

с

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

Рис. 4. Влияние установившегося вращения на положение точки отрыва на профиле

различных угловых скоростях вращения относительно х0 = 0,25 . Видно, что значения коэффициента подъемной силы совпадают очень хорошо. Зависимости для момента тангажа показаны в более крупном масштабе, в котором заметно некоторое количественное расхождение результатов. Но качественно они близки. Нелинейные изменения в зависимостях тг на относительно малых углах атаки, полученным по программе Хбэй, связаны с тем,что в этой методике такое влияние оказывает пузырь в пограничном слое на передней кромке профиля, развивающийся при переходе от ламинарного к турбулентному течению.

Расслоение аэродинамических

характеристик на малых углах атаки в зависимости от величины О хорошо описывается формулами линейной теории (3). Для центровки Хо = 0,5 зависимости т2 (а), полученные в расчетах при различных О, совпадают между собой на малых а при отсутствии отрыва потока. При центровке Хо = 0,75 аналогичная картина наблюдается для зависимостей су (а) . Результаты численных расчетов для больших

углов атаки показывают, что вращение с положительной угловой скоростью относительно Х0 = 0,25 смещает начало отрыва потока к передней кромке профиля. При вращении относительно задней центровки х0 = 0,75 картина обратная — такое же вращение приводит к смещению координаты точки отрыва назад.

По результатам расчетов с помощью программы ХРэй была вычислена функция влияния угловой скорости вращения профиля на положение точки отрыва на его верхней поверхности йх5/йО в зависимости от угла атаки для различных центровок. Эти результаты показаны на рис. 4. С помощью данных результатов и формул (4), (5) можно оценить зависимость вращательных производных аэродинамического демпфирования с^° и т° от угла

атаки. Результаты этих расчетов представлены на рис. 5. Там же показаны результаты вычисления этих производных по данным расчетов с помощью программ ЦАГИ и ХРэй. Видно, что при развитии отрыва потока производные демпфирования изменяются в зависимости от угла атаки нелинейно. В целом оценки по линеаризованной теории тонкого профиля с учетом отрыва потока позволяют описать наблюдаемые в расчетах зависимости при различных положениях центровки в широком диапазоне углов атаки.

4. Таким образом, проведенные численные исследования показывают, что расчеты аэродинамических характеристик профиля в широком диапазоне углов атаки, включающем отрывные режимы обтекания, по программам ЦАГИ и Хбэй дают результаты, очень близкие между собой, хотя в программах использованы различные методики, в том числе и по переходу от ламинарного течения к турбулентному. Это справедливо как для прямолинейного движения профиля, так и для случаев его установившегося вращения по тангажу.

Наличие вращения профиля приводит к изменению его аэродинамических

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

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

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

В статье не приводится сравнение расчетных результатов с данными эксперимента, потому что для получения последних необходимы специальные динамические установки — «ротативные машины». Авторам неизвестны опубликованные результаты для двумерных профилей, полученные на таких установках. Существующие же в ЦАГИ динамические установки по вынужденным колебаниям моделей в аэродинамических трубах дают не просто производные аэродинамических коэффициентов по углам скорости тангажа Q, а их комбинации с производными по параметру а .

ЛИТЕРАТУРА

1. Храбров А. Н. Неединственность ламинарного отрывного обтекания профиля под углом атаки в схеме Кирхгофа // Ученые записки ЦАГИ.— 1985. Т. XVI, № 5.

2. Колинько К. А., Храбров А. Н. Математическое моделирование нестационарной подъемной силы крыла большого удлинения в условиях срыва потока //

Ученые

записки ЦАГИ.— 1998. Т. XXIX, № 3 — 4.

3. Седов Л. И. Плоские задачи гидродинамики и аэродинамики / Изд. третье, переработанное.— М.: Наука.— 1980.

4. Wolkov A. V., Lyapunov S. V. Application of viscous-inviscid interaction methods for a separated flow calculation about airfoil and high-lift systems // ICAS Proceedings,

ICAS.— 1996.

5. Drela M. XFOIL: an analysis and design system for low Reynolds number airfoils //

Conference on Low Reynolds Number Airfoil Aerodynamics.— University of Notre Dame.—

June, 1989.

6. Drela M., Giles M. B. Viscous-inviscid analysis of transonic and low Reynolds number airfoils//AIAA J.— 1987. Vol. 25, N 10.

Рис. 5. Численный расчет вращательных производных c? e m? и их оценка по линеаризованной теории

7. Гуржиенко Г. А. Метод искривленных моделей и его применение к изучению криволинейного полета воздушных кораблей // Труды ЦАГИ.— 1934. Вып. 182.

8. Крамер В. В. О методе искривленных моделей / Сб. «Аэродинамика неустановившихся движений».— Труды ЦАГИ.— 1960. Вып. 775.

9. Федорова И. Б. Нелинейные компоненты нормальной силы и продольного момента тел весьма малого удлинения при установившемся вращении по кругу // Труды ЦАГИ.— 1964. Вып. 940.

Рукопись поступила 25/Х 2002 г.

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