Научная статья на тему 'Новые возможности крутильно-колебательного метода Швидковского Е. Г. : идентификация реологической принадлежности среды'

Новые возможности крутильно-колебательного метода Швидковского Е. Г. : идентификация реологической принадлежности среды Текст научной статьи по специальности «Физика»

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

Аннотация научной статьи по физике, автор научной работы — И В. Елюхина, Г П. Вяткин, В П. Бескачко

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

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

Текст научной работы на тему «Новые возможности крутильно-колебательного метода Швидковского Е. Г. : идентификация реологической принадлежности среды»

УДК 532.137.3

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

И.В. Елюхина, ГЛ. Вяткин, В.П. Бескачко

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

Введение

Крутшъно-колебателъный метод [1] широко используется в вискозиметрической практике с момента его обоснования Швидковским Е.Г для ньютоновской жидкости до настоящего времени и получил значительное усовершенствование в аппаратном отношении. Теоретическая база метода позднее была развита в [2, 3] на случай неоднородных жидкостей и учета магнитогидроди-намических эффектов, ограничиваясь, однако, по-прежнему ньютоновским поведением среды. Возможности метода в изучении неньютоновских жидкостей практически не исследованы. Единственным исключением служит работа [4], где обсуждается измерение вязкоупругих свойств, причем главное внимание уделяется режиму вынужденных высокочастотных колебаний, тогда как в вискозиметре Швидковского реализуется случай свободных затухающих колебаний

низкой частоты КГ1...!О1 Гц).

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

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

Заметим, наконец, что ошибочная интерпретация реологической принадлежности среды может приводить к противоречиям в экспериментальных данных, вопрос о причинах которых, в частности, в металлических расплавах, является открытым и представляет собой предмет интенсивных дискуссий (см., например, [5]).

Математическая формулировка задачи

Пусть цилиндрический сосуд с внутренним радиусом Я и моментом инерции К подвешен вдоль своей оси на упругой нити и совершает вокруг этой оси крутильные колебания с периодом г0 и декрементом затухания колебаний <50 . При заполнении цилиндра жидкостью, во-первых, вследствие увлечения ее движущимися ускоренно стенками цилиндра возрастает эффективный момент инерции подвесной системы и увеличивается период колебаний г: г > Тд, а во-вторых, растет скорость затухания колебаний, т.е. декремент затухания 8 : 8 > 8$, вследствие дополнительной диссипации механической энергии системы, обусловленной вязким трением между подвергаемыми сдвигу слоями жидкости. Влияние каждой из этих причин на движение цилиндра определяется полем скорости течения, возбуждаемого в жидкости стенками цилиндра. Задача

Елюхина И.В., Вяткин Г.П., Бескачко В.П.

Новые возможности крутильно-колебательного метода Швидковского Б.Г.; идентификация реологической...

заключается в предсказании закона колебаний а = a(t), т.е. зависимости от времени t углового смещения а цилиндра, заполненного исследуемой жидкостью с массой M,\ и является сопряженной: движение сосуда непосредственно связано с возбуждаемым им движением жидкости.

Характер течения жидкости при прочих равных условиях зависит от ее реологического типа, который определяется реологическим уравнением состояния, представляющим собой математическую формулировку предположений, касающихся механического поведения среды. Так, для несжимаемых ньютоновских сред достаточно указать два реологических параметра: кинематическую вязкость v и плотность р, а для неньютоновских к ним добавляются, например, динамическая жесткость G (вязкоупругая жидкость), предел текучести а 0 (вязкопластическая жидкость), показатель m и постоянная Kv степенного закона (псевдопластичные и дилатантные жидкости).

Вискозиметрические методики по измерению свойств ньютоновских сред и возможности наблюдения слабоупругих свойств линейных вязкоупругих жидкостей были обсуждены ранее (см., например, [1-4, 6, 7]). Вискозиметр, заполненный такими средами, после завершения переходного процесса совершает регулярные изосинхронные колебания, характеризуемые единственным набором параметров: периодом г и декрементом затухания S . Если реологическое уравнение состояния включает пластические или нелинейно-вязкие составляющие, то возможно нарушение подобного асимптотического поведения колебаний. Изучение качественных особенностей колебательных процессов, возникающих при заполнении вискозиметра вязкопластичными и нелинейно-вязкими средами, является целью настоящей работы и здесь будет проведено на примере моделей Бингама (вязкопластичные среды), Оствальда - Вейля (среды, подчиняющиеся степенному реологическому закону: псевдопластичные (0 < m < 1 ) и дилатантные (m > 1 )) и модели, представляющей комбинацию этих типов -модели Балкли - Гершеля (нелинейные вязкопластичные среды).

Точное решение рассматриваемой задачи отсутствует даже в простейшем случае ньютоновской жидкости ввиду нелинейности уравнений ее движения, пространственного и нестационарного характера течения. Известные (см., например, [1, 8]) аналитические решения найдены с использованием ряда приближений, главными из которых являются:

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

2) скольжение между жидкостью и внутренней поверхностью цилиндра отсутствует;

3) амплитуды колебаний малы, что в свою очередь позволяет допустить:

3.1) единственной существенной компонентой поля скорости является азимутальная,

3.2) течение жидкости в цилиндре осесимметричное.

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

При сделанных предположениях математическую модель вискозиметрических экспериментов по изучению неньютоновских свойств в цилиндрической системе координат можно представить в виде:

1 ) уравнение движения жидкости

3V даГ(р 2 агф ot or г

где V - азимутальная компонента скорости, аГ(р - г<р - ая компонента тензора напряжений а ;

2) уравнение движения цилиндра

2 _

д а 25с\ да -— +---+

ât ât

lr0 J

« = (2)

к.

;вели-

где Р - момент сил, приложенных к цилиндру со стороны жидкости, Р = -2Мр 1аг(р

чины Р, К и М отнесены к единице длины цилиндра;_

Серия «Математика, физика, химия», выпуск 3 <109

3) начально-краевые условия для (1,2):

da

V{ry 0) = 0, V(R,t) = R-

dt

К(0,/) = 0, a(0) = a0,

¿/a dt

0;

(3)

f=0

4) реологическое уравнение состояния, устанавливающее соответствие между тензором напряжений с и тензором скоростей деформации В :

4.1) модель Бингама [9]

сг

7 +

гср

£0 Z)

А^-о'

.2) модель Оствальда - Вейля [10]

лри сг > сг0,

при cr<(jQ,

_ _ ^ л772-1 п

и Г(р

(4)

(5)

4.3) модель Балкли -Гершеля [9]

^г (р

Г Л^"1 4.

А-® пРи ^^О

- —(6)

Dr(p - 0 при а < ctq .

В (4)~(6) rj~vp~ динамическая вязкость, Dr(p - гср- ая компонента D, сг и Z) - вторые инва-

рианты а и D, D

D

г(р

, 0" =

¿К V

. Заметим, что для ньютоновской среды

(7

г(р

Фг<р•

(7)

Метод решения Бингамовские жидкости

Численно трудно промоделировать поведение «идеальной» Бингамовской модели и на практике удобно использовать модель где движение нетекучей компоненты рассматривается как движение ньютоновской среды, вязкость которой намного превосходит вязкость текучей компоненты:

сг

гср

vpDr(p, v

VQ+3_ ПрИ

Dp p(yr - v0)

vr = krVQ при D <

(8)

Piyr -v0)

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

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

После введения безразмерных переменных

U =

V

Т = q0t

А =

MR

=R/d, 4 = r!d, Вт

_ ст0

D = D/qQ, (9)

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

ёцц 2 К г'0/^0

где с1 = д/уо/до - толщина пограничного слоя, дд = 2п / г о - циклическая частота колебаний

пустого тигля, Вт - число Бингама, А - отношение моментов инерции «замороженной» жидкости в тигле и пустой подвесной системы, система уравнений (1-3, 8) имеет следующий вид: при Т) < ВтЦкг-\)

ди

дТ

d2U 1 dU U

(10)

Елюхина И.В., Вяткин ГЛ., Бескачко В. П.

Новые возможности крутильно-колебательного метода Швидковского Е.Г.: идентификация реологической...

д2а да 5,

+

дТ2 дТ я

+ а =

4 Ак.

дЦ Ц д{ £

при И > Впг!{кг -1)

ди

дТ

( *г

д^Ц 1 ди Ц

2 Вт . +-эщп

'ди ил

д2а да 4А +--- + а - —-

г ди Vл

дТ2 дТ л

Я

Мо

ААВт . -З-зщп 90

'ди С/^

и*

Щ4,0) = 0, Щ4о,Г) = , ЩО,Т) = О, а(0) = «о,

и± аТ

= 0.

(11)

(12)

(13)

(14)

7=0

Для решения системы дифференциальных уравнений (10)—(14) используем метод прямых. На промежутке хе[0,<^0] выберем эквидистантную сетку узловых точек: х] = у'й, к =

у = 1,2 ..(л-1), = 0, Обозначим приближенные решения в этих точках как

и} ~1](х],Т) и аппроксимируем производные по координате х трехточечными разностными формулами. Тогда при В} < Вт!(кг -1) для системы (10), (11), (14) получаем:

йи

Л

с1и

(и2 -2щ +м0) + Х7— ("2 -ио)--%

И2 ^ ' ' 2Ихх ™ -2

, "о =0;

ж

с1и

л-1

Л

1 и

(ик~2ип_х+ип_2)+—-(ч ~ип-2)—Т1

сЛ/

П 1Пхп-\ хп_1

л

Ш К ¿¡£

ГЧ~ип-\ икЛ

к

к

где

¿ = 0: ип= а0=6° = 0,1047 рад, и] = 0, (у = 1,2...(и-2),(и-1),(и + 1))

(15)

, у = 1,2...(и-1).

Для системы (12)—(14) (при >Вт/(кг-1)) преобразования проводятся аналогично, например, последнее уравнение в (15) будет иметь вид

<4+1 _ Л

_ <5о п_ 4Л - г/л--и„+1, - -—г-

/ \ ик - ^

+ Вт ■ sign

к

(16)

Полученная система обыкновенных дифференциальных уравнений (15), (16) интегрировалась методом Рунге - Кутта четвертого порядка с контролем точности и автоматическим выбором шага по времени.

Нелинейно-вязкие жидкости

Путем введения безразмерных параметров (9) и положив, кроме того, Ь - сц™ 1, где с — Ку /(уд/?), - номинальная вязкость, систему (1)—(3), (5) можно представить как

dU - Ъ dU и m—1 2 'dU сЛ + т ' d2U 1 9U и) ---+ __

дТ 4 £ [де d W

да да So

-- +--± + а

дт дТ к

4 АЬ ' ъи и] dU и т-1

& 84 £

(17)

МО

с начально-краевыми условиями (14).

Численная реализация процесса решения проводится аналогично описанному для вязкопла-стичной среды. При расчетах с т<\ необходимо вводить адаптивную неравномерную сетку вблизи оси цилиндра или использовать аппроксимацию

у = < при £>>% (18)

при П<О0,

где - пороговое значение В9 которое можно принять ~ 10 6...10 5. Результаты и обсуждение

В дальнейшем под периодом колебаний будем понимать величину г = 2ДГг, где АТТ - разность между двумя соседними моментами времени, когда угловое смещение а обращается в

нуль, а декремент затухания определим как 8 = 2 In

а\.

а2

, где а\9а2 - соседние экстремальные

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

Используемые ниже для сравнения результаты, отвечающие ньютоновской жидкости, рассчитаны по моделям (10)—(14) и (17) при т = 1 и Вт = 0 и соответствуют результатам, полученным с помощью точного решения, приведенного в [1].

Бингамовская жидкость

На рис. 1 показано поведение периода и декремента затухания колебаний как функций номера N колебания в зависимости от различных чисел Бингама Вт . Здесь тнъют и 8ньют - период и декремент затухания при Вт =Ю, а ттв - период колебаний вискозиметра, полностью заполненного твердым ядром. По мере затухания колебаний, т.е. роста N, область, где сдвиговые напряжения превосходят предел текучести oq , будет распространяться от оси цилиндра к его стенкам. Когда твердое ядро заполнит весь объем, эффективный момент инерции системы достигнет своего наибольшего значения, равного сумме моментов инерции «замороженной» жидкости Кж и К, а

вместе с ним наибольшим станет и период, так как г ~ (К + Кж)^2 . В этом состоянии 8 должен стать минимальным и совпадающим с ¿>q ввиду отсутствия диссипации механической энергии вследствие вязкого трения. Эти качественные соображения целиком подтверждаются в расчетах.

Зависимость г = t(N) (рис. 1а) является монотонно возрастающей и тем быстрее достигает

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

значения ттв9 чем больше Вт , а 8 = 8(N) (рис. 16) сначала растет до некоторого значения 8т, а затем сравнительно быстро падает до 8 ~8q. При уменьшении числа Вт кривая 8 = 8(N) смещается к оси ординат и становится монотонно убывающей. Заметим, что при малых Вт, когда колебания затухают быстрее, чем достигается максимум 8т> поведение 8(N) и t(N)/tq сходно с характерным для псевдопластичных сред (см. рис. 3 и пояснения к нему).

Число колебаний Nme, которое успевает совершиться прежде, чем пластическое ядро заполнит весь объем, можно определить как 8(Nme) = Sq или t(Nme) s ттв . При фиксированном

число N

те

зависит от находящегося в распоряжении экспериментатора отношения

А = МЯ /(2К) (рис. 2 для = 12, Вт = 1) и наоборот, что позволяет сделать приближенную

Елюхина И.В., Вяткин Г.П., Бескачко В.П.

Новые возможности крутильно-колебательного метода _Швидковского Е.Г.: идентификация реологической...

г

1.08 1.06 1.04 1.02

1-Вт «2.0, 2 - Вт = 1.0, 0.4,

4- В>п - 0.2

нъют

Л-0.2, £0-12

10

15

N

а)

0.160

0.110

0.052

0.000

.Л= 0.2, «12 1-Я« = 2.0,

\ 2-ЙЯ-1Д

Х^^ \....... - \-------3 - Вт = 0.4,

1 \ \ \ 1 \ \ \ { \ \ \ 1 \ \ 4-5га = 0.2 \ 1 \ !

\ \ 1 \ \ " \ 1.........~ V 1 1 1 \ \ \ ' \ \ \ \ 1 \ Р \ | , \ [ \ V

\ \ \ \

\1; \2 -1—!-1- \з I \4

10

б)

15

N

Рис. 1. Зависимость периода (а) и декремента (б) затухания колебаний

оценку значения числа Вт . Число Ытв растет с уменьшением А, период тнъют и декремент 5т падают с ростом и уменьшением А, и соответствующим образом изменяются тте/тньют (рис. 2) и (<Ут-<У0). Это дает возможность изменять условия наблюдаемости и идентифицируемости пластических свойств путем варьирования временного интервала изменения параметров колебаний и их чувствительности к числу Бингама.

0.05 0.10 0.15 А Рис. 2. Зависимости от А

Итв и г

те ' * ньют

Нелинейно-вязкие жидкости

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

падает, а для псевдопластических - растет. Напомним, что в ньютоновском

К

д% -у

/ дг /г

случае период колебаний т уменьшается при уменьшении вязкости у , т.е. росте ¿¡о при фиксированных Я и до > а поведение декремента 8 зависит от величины : при > с Рос"

том <5л декремент уменьшается, а при < ~ растет, где = 4,2 [1]. Эти законо-

'шах

мерности качественно подтверждаются и в случае нелинейно-вязких сред, что проиллюстрировано на рис. 3, 4. Здесь

■1

ЫУ

л

в

т-1

¿0> 0 =

аи/эг~%

(19)

где Ь определено при амплитудном значении а, а отвечает номинальной величине вязкости. В общем случае О находится как интегральная характеристика по полупериоду, и ввиду того, что четко уловить границу £од<у\я Для нелинейно-вязких сред нельзя, она колеблется в не-

больших пределах в зависимости от А и т . Отметим также, что для ньютоновских сред при больших значениях ^ с ростом период уменьшается слабее, чем декремент, что для нелинейно-вязких сред можно проследить на рис. 3: выход на асимптотический режим, характери-

зуемый неизменными во времени значениями г и 8, происходит быстрее по периоду (рис. За), чем по декременту (рис. 36). Для жидкостей со степенным реологическим законом характерны высокие значения > и в этом случае период в процессе колебаний практически не изменяется. Рис. 3, 4 построены для т >1, а при т <1 зависимости периода и декремента отображаются практически симметрично относительно прямой, соответствующей ньютоновской среде с т = 1 (зависимость 8 = 8{Щ с учетом ¿¡о*гу).

Рис. 4. Зависимость 8 - ¿>(£п »л/ )

Рис. 3. Зависимость периода (а) и декремента (б) затухания

колебаний от их номера для сред Оствальда-Вейля 1) ^ = б , т = 2; 2) = 8, т = 1.5

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

Заключение

Таким образом, в настоящей работе:

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

1.1) колебания становятся неизосинхронными - период является функцией номера колебаний до момента, пока не исчезнет текучая компонента, в случае бингамовских сред или асимптотически стремится к некоторому пределу для нелинейно-вязких сред (напомним, что последняя ситуация возможна и для вязкопластичных сред при малых числах Бингама),

1.2) декремент затухания также изменяется в процессе колебаний на всем временном интервале, где среда обладает свойством текучести, и может иметь схожий характер для бингамовских и нелинейно-вязких сред;

2) полученные результаты демонстрируют возможность экспериментальной идентификации реологической принадлежности жидкости:

8 И N

Рис. 5. Зависимость периода и декремента затухания колебаний от их номера для моделей: Бингама---, Оствальда-Вейля....... , Балкли-Гершеля

А = 0.2, = 12, т = 2, Вт = 0.4

Елюхина И.В., Вяткин Г./7., Бескачко В.П.

Новые возможности крутильно-колебательного метода Швидковского Е.Г.: идентификация реологической...

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

2.2) для нелинейно-вязких жидкостей со степенным реологическим законом режим колебаний с единственными значениями т и S устанавливается не за конечное число колебаний, а лишь асимптотически, а их принадлежность к дшатпантным или псевдопластическим средам можно выявить по характеру зависимостей г = t(N) и 8 = 8{N).

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

Работа выполнена при поддержке Минобразования РФ и Правительства Челябинской области (грант 2003 г. для молодых ученых вузов Челябинской области).

Литература

1. Швидковский Е.Г. Некоторые вопросы вязкости расплавленных металлов. - М.: ГИТТЛ, 1955. - 206 с.

2. Бескачко В.П., Вяткин Г.П., Писарев Н.М., Хисматулин М.Б. Теория крутильного вискозиметра, помещенного в осевое магнитное поле // Магнитная гидродинамика. - 1992. - № 2. -С. 65-70.

3. Бескачко В.П., Хисматулин М.Б., Щека А.И. Теория крутильного вискозиметра, заполненного стратифицированной жидкостью // Физико-химические основы металлургических процессов. - Челябинск: ЧГТУ, 1992. - С. 47-52.

4. Kleiman R.N. Analysis of the oscillating-cup viscometer for the measurement of viscoelastic properties//Phys. Rev., 1987.-Vol. 35, № l.-P. 261-275.

5. Островский О.И., Григорян B.A. О структурных превращениях в металлических расплавах // Изв. вузов. Черная металлургия. - 1985. - № 5. - С. 1-12.

6. Елюхина И.В. Планирование оптимального эксперимента по одновременному определению вязкости и плотности ньютоновской среды // Вестник ЮУрГУ. Серия «Математика, физика, химия». Вып. 1.-Челябинск, 2001.-№ 7 (07) - С.71-76.

7. Елюхина И.В. К вопросу наблюдаемости упругих свойств жидких сред в вискозиметриче-ском эксперименте по Швидковскому Е.Г. // Там же - С.77-81.

8.Kestin J., Newell G.F. Theory of oscillating type viscometers: the oscillating cup. Part I. // J. Appl. Math. Phys., ZAMP, 1957. - V. 8. - P.433^149.

9. Balmforth N.J., Craster R.V. A consistent thin-layer theory for Bingham plastics // J. Non-Newtonian Fluid Mech., 1999. -№ 84. -P.65-81.

10. Fang P., Manglik R.M., Jog M.A. Characteristics of laminar viscous shear-thinning fluid flows in eccentric annular channels // J. Non-Newtonian Fluid Mech., 1999. - № 84. - P. 1-17.

Поступила с редакцию 15 апреля 2003 г.

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