Научная статья на тему 'Эволюция профилей скорости и турбулентной вязкости в системе течений со сгонно-нагонным и плотностным потоками'

Эволюция профилей скорости и турбулентной вязкости в системе течений со сгонно-нагонным и плотностным потоками Текст научной статьи по специальности «Физика»

CC BY
44
8
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
СТРАТИФИЦИРОВАННЫЕ ТЕЧЕНИЯ / STRATIFIED CURRENTS / ТУРБУЛЕНТНЫЙ ОБМЕН / TURBULENT EXCHANGE / СГОННО-НАГОННЫЙ ПОТОК / WIND-INDUCED FLOW / ГИДРОДИНАМИЧЕСКАЯ УСТОЙЧИВОСТЬ / HYDRODYNAMIC STABILITY / ТУРБУЛЕНТНАЯ ВЯЗКОСТЬ / TURBULENT VISCOSITY

Аннотация научной статьи по физике, автор научной работы — Самолюбов Борис Исаевич, Иванова Ирина Николаевна

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

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

Текст научной работы на тему «Эволюция профилей скорости и турбулентной вязкости в системе течений со сгонно-нагонным и плотностным потоками»

60

ВМУ. Серия 3. ФИЗИКА. АСТРОНОМИЯ. 2014. № 5

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

Б. И. Самолюбовa, И. Н. Ивановаb

Московский государственный университет имени М. В. Ломоносова, физический факультет, кафедра физики моря и вод суши. Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2.

E-mail: a [email protected], b [email protected] Статья поступила 02.06.2014, подписана в печать 16.06.2014.

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

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

УДК: 556.532.5. PACS: 92.10.Lq, 92.10.Wa, 92.40.-t, 92.40.Ni.

Введение

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

1. Объект и методика исследований

Анализируемые результаты получены 19.09.2007 г. в Петрозаводской губе Онежского озера в период слабой стратификации перед осенним этапом полного выравнивания температур воды по всей глубине [4]. Вертикаль зондирований с глубиной 23 м располагалась в центре восточной части залива. Средняя по глубине частота Вяйсяля-Брента N составляла (2-3) • 10~2 Гц. Регистрации профилей скорости течения по глубине водоема с дискретностью 10 мин при шаге по вертикали 0,5 м велись доплеровским профилографом Я0СР-600 (Лапёегаа). Прибор ЯЭСР-600 в режиме донной постановки обеспечивал также регистрацию уровня свободной поверхности. Одновременно зондом ЯСМ 9 LW (той же фирмы) измерялись профили скорости и, температуры Т и электропроводности Ск, концентраций кислорода О2 и мутности Ти. Точности измерений

а

0.0001 0.0002 Р, г/смз

А

-60

-20

0

z, м

Z, м

¿0 _Z'M 20

а 16

12q Si 12

8 О - тУ 8

4о ^О 4

о

U, см/с

0

20

40

60

20 16 12 8 4 0

Ки, см2/с ^

7 S, мг/л

Рис. 1. Распределения по всей глубине: а — скорости течения и (г) и изменения плотности воды с глубиной Др(г) по данным измерений (точки, пунктир) и расчета и по моделям течения (сплошная линия); б — коэффициента турбулентного обмена; в — измеренной концентрации взвеси Ти (точки) и рассчитанной

по модели (сплошная линия)

и, Т, Ск, 02 и Ти составляли 0.5 см/с, 0.02°С, 0.02 мСм/см, 0.25 мг/л и 0.4 ЫТи. Параллельно анеморумбометром М-63М измерялась скорость ветра с точностью ± 0.25 м/с.

2. Структура течения

Представление о структуре течения дают профили и распределения по глубине и во времени величин и, Т, 5р, 5, а также коэффициента обмена Ки, рассчитанного по представленной ниже методике (рис. 1, 2). Здесь 5 — концентрация взвеси, находившаяся по результатам калибровки турбидиметра на основании весового анализа проб воды, 5р = р(г) — р(Н — 0.5 м) — изменение плотности с глубиной, Н — глубина места. В серии зондирований на срочной станции было получено 34 профиля и , сформированных под влиянием сгонно-нагонного и плотностного течений [5]. Сгон-но-нагонное течение содержало дрейфовый (вверх по заливу) и компенсационный потоки. Средние скорости составляли 30 см/с для дрейфового течения и 6 и 8 см/с для компенсационного и плотностного потоков, направленных вниз по заливу в открытое озеро (положительные значения и (г) на рис 1, а). Распределение 5р определялось главным образом термической стратификацией при разности температур воды у поверхности и у дна порядка 2.5° С. Изменение 5р с глубиной характеризуется максимумами вертикального градиента не только в термоклине, но и в придонном слое толщиной порядка 1 м. Профиль турбулентной вязкости Ки содержит пики в дрейфовом потоке и на верхних границах компенсационного и плотностного течений. Распределение концентрации взвеси 5 в целом устойчивое (рис. 1, в).

В поле скорости (рис. 2, а) колебания изотах определяются влиянием внутренних волн, которым в нижнем 10-м слое соответствуют периоды и высоты 1-3 ч и 2-5 м. С аналогичными интервалами времени происходят вспышки интенсивности турбулентного обмена, которые характеризуются, отмеченными выше пиками Ки (рис. 2, б). В поле температуры прослеживаются вызванное ветром заглубление термоклина и волновые колебания изотерм в его области.

3. Распределение скорости в системе течений

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

U (z) = Uc + Udc + Ued.

(1)

Вертикальное распределение скорости для компоненты Udc получается на базе решений, полученных в работах [6, 7]. Плотность воды р рассчитывается по данным измерений параметров состояния воды по уравнению Чена-Миллеро [8].

Распределение Uc получалось из уравнения Рей-нольдса для сгонно-нагонного потока над термоклином (как над жидким грунтом), высота которого над уровнем дна z^ = 1 -6 м совпадала с верхней границей придонного плотностного течения (рис. 1,а). Уровню zig соответствовала изотерма 10° С(рис. 2, в). Оценки придонного и приповерхностного экмановских слоев,

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

drXz d(s

ИГ = -gpa1 dX'

(2)

в котором тхг = рКи(г)Щц — турбулентное напряжение, Сх — возвышение уровня воды относительно равновесного положения. Распределение турбулентной вязкости Ки представлено в следующем разделе. Здесь а1 — коэффициент, учитывающий влияние сил, не включенных в уравнение динамики, в том числе обусловленных воздействием крупномасштабных внутренних волн и трехмерностью структуры течения. Среднее значение а1 составляет 0.83 (± 10%) по результатам проверки модели.

Уравнение (2) интегрируется с граничным условием Тхг\г=Сх = тw, где тш = а0расодиш\иш\ — напряжение трения ветра, С0А = 1.4 • 10~3 — коэффициент сопротивления на границе вода-воздух (при регистрации скорости ветра и^ в 10 м над поверхностью воды), ра = 1.2 • 10-3 г/см3 — плотность воздуха при Т = 20 ° С, аа = 2.9 — коэффициент, зависящий от степени развития ветрового волнения [10].

Интегрирование (2) по высоте г в слое г^ < г < Н дает ис в виде

Uc = ^ Ро

dz K

a1g

дС

dx

H + Zs - (z + zig)

Ku

dz. (3)

zig

zig

Возвышение уровня воды регистрировалось дат-

чиком давления профилографа ЯЭСР 600, работавшего в режиме донной постановки.

Для вихревой составляющей Ц^ принято за основу теоретическое распределение для возмущений скорости, вызванных развитием завихренности в пограничном слое [11]. Изменения такой природы на профиле скорости отмечены вблизи уровня смены знака и(г) при переходе от дрейфового потока к компенсационному течению. Согласно результатам, представленным в [12], подобные вихревые структуры, порожденные развитием неустойчивости, типичны для аналогичных зон перехода скорости через нулевое значение. Причем предполагается, что слой нулевой скорости ведет себя подобно слою разрыва в поле плотности, вблизи которого генерируется завихренность. Этот факт подтверждается в [12] наблюдениями аналогичных вихрей в атмосфере. По [11] профиль ие^ характеризуется квазигармонической функцией с амплитудой, экспоненциально затухающей с удалением от зоны генерации завихренности. Центр этой зоны располагается в точке перегиба при г = ге^ на профиле и (г). На этой высоте течение теряет устойчивость. Неустойчивость вызывает генерацию завихренности, горизонтальная компонента скорости которой записывается по [11] в виде

Ued = Are-kelz-zedl cos(kr\z - zed\ + фг),

(4)

ВМУ. Серия 3. ФИЗИКА. АСТРОНОМИЯ. 2014. № 5 а

Н, м

13:30 16:00 - 18:30 ч

б

Ки, см/с

14 16 18 1,ч

Рис. 2. Распределения скорости течения (а), коэффициента турбулентного обмена (б) и температуры воды (в

по высоте над уровнем дна и во времени

где Аг, кг, фг — максимальная амплитуда, волновое число и фаза, ке — коэффициент затухания. Волновое число кг = 2п/(А2)еА, где (А2)^ — вертикальная длина волны для возмущений, равная удвоенной глубине центра слоя нулевых скоростей. То есть в пределах этой глубины укладывается половина длины волны. Коэффициент ке оценивался в соответствии с (4) как значение \г — 2е^\-1, при котором амплитуда возмущений убывала в е раз относительно ее значения на уровне смены знака скорости течения. Для максимального значения амплитуды Аг получено выражение Аг = 0.34-(Н/гтоМ1 — 2(гт/Н))• (Ц^о. Значения перечисленных параметров составляли: кг = кг0 = 0.57 м -1, фг = 0.87п, к-1 = 0.15 Н, 2Т — высота термоклина над уровнем дна, Ц5ЦГ( — приповерхностная скорость течения, которая более подробно рассматривается ниже. Подстрочный индекс ноль соответствует значениям параметров при I = 0.

Профиль ит под уровнем жидкого грунта 2^ в слое смешения плотностного потока (2 = 2т ч—2^) определяется как ит = 1\д2Ц\, где I = 0.06 • (2^ — гт), \дги\ =

= (6Цт/(% — 2т))^( 1 — £) , Цт = Цт (Ар, 2^, ¿з, С0) [7].

Компоненты ЦС и Ц]С гладко сшиваются при 2 = 2^, а с ЦС в слое нулевых скоростей (в окрестности точки перегиба при 2 = 2е^ на профиле Ц(2)). Приведенная методика проверена в расчетах профилей скорости по (1), (3), (4) и концентрации взвешенных частиц по уравнению диффузии взвеси [13] (рис. 1).

4. Распределение турбулентной вязкости и его эволюция

В формуле (3) применялось выражение коэффициента турбулентного обмена из [13] в виде

Ku = U +

VeC) i/ (-

0.4 N 2l2 \

u2 + efc )

1/2

(5)

где uT = l\dzU\ — сдвиговая скорость; e¡c = 4-10-4-UFD ■ ■ (1 - \dzU/(dzU)max\) — оценка энергии вертикального турбулентного обмена в слоях воды с градиентами скорости dzU, близкими к нулевым. Здесь UFD — скорость, средняя по всей глубине; N2 = (g/p)(-dzp) — частота плавучести; р — плотность воды. Профиль масштаба турбулентности (пути смешения) I представляется в виде

I =

/«V 0

\к(Н -

+ (кг/La)) при 2 < 2П

- г)/{ 1 + (к(Н - 2)/Las)) при 2 ^ 2п

(6)

где 2mc — уровень максимума скорости компенсационного течения. Здесь k = 0.4, La и Las — интегральные масштабы турбулентности для придонной (2 < 2mc) и приповерхностной (2 ^ 2mc) областей [10, 12]. Из условия сшивания распределений La и Las при 2 = 2mc следует, что

Las = k ' 2mc • (H 2mc) •La / ((H 2mc) • (La + k'2mc) 2mc 'La) .

Значение l приведено выше. Из этого условия по (6) находится выражение La в виде

La = K2mc/((k2mc/0.06 • (% - 2m)) - 1) .

Проверка правильности определения Ku проводилась путем сравнения измеренных и теоретических распре-

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

делений концентрации взвеси, полученных по методике, приведенной в [7] с применением Ки из (5).

С применением разработанной и проверенной методики расчета профилей коэффициента турбулентного обмена были получены распределения турбулентной вязкости по глубине и во времени в сериях зондирований в Петрозаводской губе Онежского озера. Диапазон значений коэффициента обмена составлял 0.1-50 см2/с. Эти результаты анализировались совместно с распределениями параметров течения и состава воды. Анализ распределения коэффициента обмена по глубине и во времени в системе со сгонно-нагонным и плотностным потоками дал следующие результаты.

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

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

1. Волновые колебания коэффициента обмена в придонном слое, близкие по периоду к изменениям толщины этого слоя, под влиянием внутренней волны. Период колебаний (со средним значением 1 ч) аналогичен зарегистрированному в полях Ц, Т и 5.

2. Увеличение (в 1.5-2 раза) максимумов коэффициента обмена в дрейфовом потоке в моменты максимальных ускорений течения. Колебания границ зон этих максимумов на распределении коэффициента обмена по глубине и во времени также имеют волновой характер с указанным выше периодом. Причем нижние границы колеблются в противофазе с термоклином. Мощнейшим из максимумов коэффициента обмена (при I = 15, 15.4, 18, 18.4 ч) в дрейфовом потоке соответствует повышение турбулентной вязкости, которое заметно не только в этом потоке, но и в придонной области. Такое проникновение турбулентности в глубинные слои сопровождается активным перемешиванием вод, которое приводит (при I = 15 ч) к локальному по времени охлаждению приповерхностных слоев.

5. Оценка распределения турбулентной вязкости при неизвестном профиле скорости

Если профиль Ц(2) не задан, то функции е^ и ит, входящие в выражение Ки (5), определяются по найденным и проверенным авторами полуэмпирическим выражениям. Эти выражения получены с учетом выявленных особенностей распределений еС и ит.

Для величины е^ проверена возможность применения приближения квазиинвариантности по глубине. При этом функция е^(2) заменяется ее средним по вертикали значением еС. Анализ зависимости еС от устойчивости стратификации с учетом пропорцио-

64

ВМУ. Серия 3. ФИЗИКА. АСТРОНОМИЯ. 2014. № 5

нальности скоростей, средних по глубине UFD и приповерхностных Usurf = U\z=H-0.5 м позволил получить выражение вида

f = 0.017 • Us2urf (102 (Npo/Nq^ — 0 .

(7)

Здесь N^0 и Nqh — частоты плавучести, средние по всей глубине и над термоклином на профиле плотности. Приповерхностная скорость течения и5ЦГ| определяется при известном возвышении уровня с применением выражения из [14] как

Usurf = kw • VgChb/CadL in (h—bHl) '

(8)

где hb — разность глубин на границах бассейна длиной L, HLA — средняя глубина, Cad = 0.0026 — коэффициент сопротивления на границе вода-воздух, kw — безразмерный коэффициент ветрового возвышения уровня. Коэффициент kw определялся по данным измерений Usurf и (s на вертикали зондирований и составлял 0.034.

Профиль uT под уровнем жидкого грунта zig в слое смешения плотностного потока ( z = zm -—zig ) определяется как uT = 0.36 • Um£(1 - £), а при z < zm профиль uT квазилинеен (см. пояснения к (6)). Над уровнем жидкого грунта при zig < z < z* профиль uT принимается квазиоднородным исходя из формы распределения uT, найденного с применением теоретического профиля U. Из-за сглаживания неоднородностей профиля uT при интегрировании (в ходе решения уравнения (3)) допустима аппроксимация uT = const при zig < z < z*, где z* = 17 м (рис. 3).

При z ^ z* на профилях uTgстойчиво воспроизводятся два максимума, которые наилучшим образом аппроксимируются распределениями, аналогичными полученным для слоя смешения в работе [7], в виде

Здесь (u^i — максимум ^ на уровнях нижнего и верхнего пиков i = 1 и 2, = \(z — z^)^/^)up —

( \ 02

(zтi)low)\, ut1 = 0.0005 • Uw • (T^vt) , ut2 = 2ut1 •

Nav \

)lo

v (Nav)0 ,

(zтi)up — высоты ниж-

— (uтm)i • ^тi(1 .

(9)

У^уЬ/ ' у ти 1о* у ти иР

ней и верхней границ для каждого из пиков I = 1 и 2. Эти границы определяются как (гт1)[0№ = г*, (гт1 )ир = (гт2)^ = 1.5 • Н • и-0 3, (гт2)ир = Н.

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

Заключение

1. Выявлены закономерности преобразований структуры системы течений со сгонно-нагонным и плотност-ным потоками.

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

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

4. Обнаружены периодически воспроизводящиеся элементы структуры распределения коэффициента обмена: а) резкие увеличения турбулентной вязкости с пиками в дрейфовом потоке; б) сопутствующие им

Z, M

U, см/с

-50 -40 -30 -20 -10 0 10

з их, см/С

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

u

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

Работа выполнена при финансовой поддержке РФФИ (грант 14-05-00822а).

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

1. Umlauf L., Lemmin U. // Limnol. Oceanogr. 2005. 50, N 5.

P. 1601.

2. De Cesare G., Boillat J.L. // XXX IAHR Congress. Auth.

Thessaloniki, Greece. 24-29 August 200З. V. 1. Theme C.

P. 381.

3. Serra T., Vidal J., Casamitjana X. // Limnol. Oceanogr.

2007. 52, N 2. P. 620.

4. Показеев К.В., Самолюбов Б.И., Филатов H.H. // Метеорология и гидрология. 2012. № 2. C. 8З.

5. Самолюбов Б.И., Иванова И.Н. // Изв. РАН. Сер. физ. 2010. 74, № 12. C. 1770.

6. Саркисян А.С. Численный анализ и прогноз морских течений. Л., 1978.

7. Самолюбов Б.И. Плотностные течения и диффузия примесей. М., 2007.

8. Chen C.T.A., Millero F.J. // Limnol. Oceanogr. 1986. 31, N 3. P. 657.

9. Horn D.A., Imberger J., Ivey G.N. // J. Fluid Mech. 2001. P. 181.

10. Зырянов В.Н., Фролов А.П. // Водные ресурсы. 2006. 33, № 1. C. 1.

11. Бетчов В., Криминале В. Вопросы гидродинамической устойчивости. М., 1971.

12. Скорер Р. Аэрогидродинамика окружающей среды. М.,

1980.

13. Самолюбов Б.И. // Вестн. Моск. ун-та. Физ. Астрон. 2012. № 4. С. 81 (Moscow University Phys. Bull. 2012. 67, N 4. P. 398).

14. McCormick M.J, Scavia D. // J. Water Resources Res.

1981. 17, N 2. P. 305.

The evolution of velocity profiles and turbulent viscosity in a system of currents with wind-induced and density flows

B.I. Samolyubova, I.N. Ivanovab

Department of Marine and Inland Water Physics, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia.

E-mail: a [email protected], b [email protected].

Abstract—The results of an in-situ study and mathematical modeling of a system of currents with wind-induced and density flows are presented. The model is based on the obtained dependences of its key parameters on the density-field characteristics, changes in the water level, and the water-reservoir topography. The revealed features of the of shear-velocity profile are indicated and a function of its distribution from the bottom to the open surface of the water is given. The model that was developed is verified by the data from measuring the parameters of currents and composition of water in the Petrozavodsk Bay of Onega Lake in September 2007. The regularities for the evolution of distributions of the coefficient of the turbulent exchange and current velocity over depth and in time are revealed.

Keywords: stratified currents, turbulent exchange, wind-induced flow, hydrodynamic stability, turbulent viscosity. PACS: 92.10.Lq, 92.10.Wa, 92.40.-t, 92.40.Ni. Received 2 June 2014.

English version: Moscow University Physics Bulletin 5(2014). Сведения об авторах

1. Самолюбов Борис Исаевич — доктор физ.-мат. наук, профессор, гл. науч. сотрудник; тел.: (495) 939-10-46, e-mail: [email protected].

2. Иванова Ирина Николаевна — канд. физ.-мат. наук, ст. науч. сотрудник; тел.: (495) 939-10-46, e-mail: [email protected].

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