Экономико-математическое моделирование: Колонна для перегона коньячного спирта, Курсовая работа

Содержание

Вступление

1  Краткое описание колонны перегонки коньячного спирта

2  Математическая модель установки и преобразование ее в пространство состояний

3  Преобразование математической модели в дискретное время и ее проверка с помощью построения разгонных характеристик.

4  Синтез многомерного ПИ-регулятора

5  Моделирование замкнутой системы и оценка качества переходных процессов

6  Преобразование модели регулятора в форму, отвечающую ее реализации в программном обеспечении

7  Выбор технических средств реализации системы управления

Выводы

Литература

Приложение1 – текст программы


Аннотация

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


Вступление

В данной курсовой работе проводится проектирование математического обеспечения колонны перегонки коньячного спирта для дальнейшего его использования в программе обеспечения АСУТП установки. Были построены характеристики объекта: кривые разгона, переходные процессы без подачи возмущений. Вышеперечисленные операции выполнялись согласно исходной матрице передаточных функций. Все математические расчеты были произведены с помощью пакета MATLAB7. Выбраны технические средства автоматизации для реализации спроектированной системы управления.


Задание

Колонна перегонки коньячного спирта

Рисунок 1.- Технологическая схема

Таблица 1. Матрица передаточных функций объекта

u1, м3/с , брага u2, кг/с, пар

y1

5+-0.1

МПа

y2,

0.7+-0.05

у1,…, у2 – управляемые переменные (измерения), u1,…,u2 – управляющие воздействия, время в секундах


1 Краткое описание колонны перегонки коньячного спирта

Колонная установка - это устройство для перегонки виноматериалов в непрерывном потоке.

Установка для получения коньячного спирта, устройство для перегонки виноматериала и коньячного спирта-сырца. Различают установки периодического и непрерывного действия (см. колонная установка). Установки периодического действия бывают однократной (ПУ-500) и двукратной (установка для получения коньячного спирта шарантского типа) перегонки. Установка ПУ-500 позволяет получать коньячный спирт крепостью 62—70% об. непосредственно из виноматериала и наряду с аппаратом шарантского типа нашла широкое распространение в коньячном производстве. Состоит из перегонного куба 1 полезной емкостью 500 дал, ректификационной колонны 2 с 3—4-колпачковыми тарелками, вертикального кожухотрубчатого дефлегматора 3, подогревателя виноматериала 4, холодильника 5, сборников дистиллята 6 и контрольно-измерительных приборов. Установка снабжена также вакуум-прерывателем, обеспечивающим безопасность работы аппарата. В перегонном кубе виноматериал доводят до кипения.

Конденсат (флегма) непрерывно подается на орошение колонны, а несконденсировавшиеся пары спирта из дефлегматора поступают в змеевик подогревателя для нагрева до 60—70°С новой порции виноматериала или непосредственно направляются в холодильник, откуда охлажденный дистиллят стекает в сборники. В процессе дистилляции отделяют головную, среднюю (коньячный спирт) и хвостовую фракции (см. перегонка виноматериалов ). Головную фракцию крепостью 80—87% об. отбирают в количестве 0,8—1,2% от безводного спирта, содержащегося в навалке, и используют для получения ректификованного спирта. Затем приступают к отбору средней фракции (коньячный спирт), крепостью 62—70% об. Процесс ведут в течение 4—4,5 ч. При снижении крепости дистиллята до 45—50°С начинают отбор хвостовой фракции, которую добавляют к виноматериалу и после пятикратного возврата выделяют и направляют на ректификацию. Дистилляцию прекращают при показании спиртомера 1% об. Кубовый остаток (барду) сливают и направляют на утилизацию. За 5—6 ч до окончания процесса в подогреватель загружают 450 дал виноматериала и 50 дал хвостовой фракции от предыдущей перегонки. Общая продолжительность процесса (загрузка куба, перегонка вина, слив остатка) 12 ч. Производительность установки (при перегонке виноматериала крепостью 10% об.) 100 дал безводного спирта в сутки. На базе ПУ-500 разработана новая установка Б2-ВУФ производительностью до 200 дал коньячного спирта в сутки.

Дефлегмация (от де… и греч. phlégma — мокрота, влага) — частичная конденсация смесей различных паров и газов с целью обогащения их низкокипящими компонентами.Дефлегмация основана на преимущественной конденсации высококипящих компонентов при их охлаждении. Дефлегмация представляет собой разновидность противоточной фракционированной конденсации. Дефлегмация пользуются как промежуточной стадией при разделении газовых смесей, а также в процессах дистилляции и ректификации. Самостоятельно дефлегмацию применяют при разделении газовых смесей, компоненты которых значительно различаются по температуре конденсации.

В странах СНГ колонные установки применяют для получения коньячного спирта с отделением эфиро-альдегидной фракции (К-5М; КПИ) и без отбора головного погона (К-5, производство Болгарии).

Основные части промышленной колонной установки:

- колонна с контактными устройствами,

- дефлегматор,

- холодильник,

- подогреватель.

В колонную установку К-5 виноматериал из подогревателя подается непрерывно на питающую тарелку, пары поступают в дефлегматор, где примерно 2/3 конденсируются, образуя флегму, стекающую в колонну. Оставшаяся часть паров конденсируется в холодильнике, а дистиллят отводится в сборник. Спирт содержит значительное количество головных примесей, снижающих его качество.

Особенностью других колонных установок являются устройства для тепловой обработки вина и для отделения эфироальдегидной фракции.

Наибольшее распространение получила колонная установка К-5М, которая (наряду с общими элементами) содержит перегреватель, где вино выдерживается 3 ч при температуре 105°—110°С, охладитель, в котором перегретое вино охлаждается до 85°—90°С, и эпюрационную колонну, обеспечивающую отбор головной фракции (1—3% в пересчете на безводный спирт). Производительность К-5М до 400 дал/сутки безводного спирта при расходе 160—180 дал/ч вина (10% об.).

2 Математическая модель установки и преобразование ее в пространство состояний

 

Математическая модель в виде матрицы передаточных функций приведена в таблице 2.

Таблица.2

u1, м3/с , брага u2, кг/с, пар

 y2,

0.7+-0.05


На рисунке 2 представлена блок – схема модели колонны.

Рисунок 2 – Блок-схема модели колонны

В исходных данных, модель дана как мы видим в виде матриц передаточных фунцый. Для преобразования передаточных функций в пространство состояний использовали соотношения. Наиболее простой аппроксимацией опоздания является замена его инерционным звеном первого порядка. Для проверки правильности преобразования следует найти собственные значения системы с помощью функции Eig и убедиться, что или все собственные значения имеют отрицательные действительные части (система постоянна), или число нулевых собственных значений совпадает с числом интегральных звеньев в исходной модели. Окончательно система должна быть представлена матрицами A,B,C,D.


 
 

Рисунок 3. Развернутая структурная схема системы с учетом запаздывания

Исходя из систем получим матрицы модели в пространстве состояний

 

где х- состояние систем;

y- измеряемые входы;

f- возмущение;

u- управление.

Матрицы системы имеют вид:

A=[-1/35 0 0 0 0 0 0 0 0 0 0 0 0 0;

0 -1/129 0 0 0 0 0 0 0 0 0 0 0 0;

0 4/48 -2/48 0 0 0 0 0 0 0 0 0 0 0;

0 0 0 -1/38 0 0 0 0 0 0 0 0 0 0;

0 0 0 4/9 -2/9 0 0 0 0 0 0 0 0 0;

0 0 0 0 0 -1/110 0 0 0 0 0 0 0 0;

0 0 0 0 0 4/134 -2/134 0 0 0 0 0 0 0;

0 0 0 0 0 0 0 -1/13.5 0 0 0 0 0 0;

0 0 0 0 0 0 0 0 -1/98 0 0 0 0 0;

0 0 0 0 0 0 0 0 4/133 -2/133 0 0 0 0;

0 0 0 0 0 0 0 0 0 0 -1/50 0 0 0;

0 0 0 0 0 0 0 0 0 0 4/12 -2/12 0 0;

0 0 0 0 0 0 0 0 0 0 0 0 -1/186 0;

0 0 0 0 0 0 0 0 0 0 0 0 4/150 -2/150];

Матрица входа:

B=[-42.5/35 0;

-72.5/129 0;

0 0;

1720/38 0;

0 0;

730/110 0;

0 0;

0 0.994/13.5;

0 0.459/98;

0 0;

0 -6.9/50;

0 0;

0 -5.1/186;

0 0];

Матрица измерений:

C=[1 -1 1 0 0 0 0 1 -1 1 0 0 0 0;

0 0 0 -1 1 -1 1 0 0 0 -1 1 -1 1];

где матрица системы:

D=[0 0;0 0];

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

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

dt=0.01/max(abs(eig(A)))

t=0:dt:999;

[Ad,Bd]=c2d(A,B,dt);

dt=0.4500

Проверить найденную модель в дискретном времени следует с помощью расчета разгонных характеристик. Для этого следует использовать функцию dstep. Для вывода графиков следует использовать функции: subplot, plot, grid.

Ad =

Columns 1 through 8

0.9872 0 0 0 0 0 0 0

0 0.9965 0 0 0 0 0 0

0 0.0371 0.9814 0 0 0 0 0

0 0 0 0.9882 0 0 0 0

0 0 0 0.1892 0.9048 0 0 0

0 0 0 0 0 0.9959 0 0

0 0 0 0 0 0.0134 0.9933 0

0 0 0 0 0 0 0 0.9672

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

Columns 9 through 14

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0.9954 0 0 0 0 0

0.0135 0.9933 0 0 0 0

0 0 0.9910 0 0 0

0 0 0.1439 0.9277 0 0

0 0 0 0 0.9976 0

0 0 0 0 0.0119 0.9940

Bd =

-0.5429 0

-0.2525 0

-0.0047 0

20.2483 0

1.9628 0

2.9803 0

0.0200 0

0 0.0326

0 0.0021

0 0.0000

0 -0.0618

0 -0.0045

0 -0.0123

0 -0.0001

Построим разгонные характеристики с помощью функций dstep, subplot, plot, grid.

Рисунок 4.Кривые разгона.

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

4 Синтез многомерного ПИ-регулятора

Для синтеза ПИ-регулятора полученные матрицы должны быть расширены в матрицы A1, B1, C1:

A1=[Ad zeros(n,l); C eye(l)];

B1=[Bd;zeros(m)];

C1=[C eye(l)];

Матрицы параметров регулятора должны быть расчитаны с помощью функции dlqr.

K=dlqr(A1,B1,Q,R)

L=dlqr(A1',C1',Q1,R1)'

Весовые матрицы Q1,R1,Q,R выбраны как единичные (для простоты матрицы генерирует функция eye).

Матрицы имеют вид:

A1 =

Columns 1 through 8

0.9872 0 0 0 0 0 0 0

0 0.9965 0 0 0 0 0 0

0 0.0371 0.9814 0 0 0 0 0

0 0 0 0.9882 0 0 0 0

0 0 0 0.1892 0.9048 0 0 0

0 0 0 0 0 0.9959 0 0

0 0 0 0 0 0.0134 0.9933 0

0 0 0 0 0 0 0 0.9672

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

1.0000 -1.0000 1.0000 0 0 0 0 1.0000

0 0 0 -1.0000 1.0000 -1.0000 1.0000 0

Columns 9 through 16

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0.9954 0 0 0 0 0 0 0

0.0135 0.9933 0 0 0 0 0 0

0 0 0.9910 0 0 0 0 0

0 0 0.1439 0.9277 0 0 0 0

0 0 0 0 0.9976 0 0 0

0 0 0 0 0.0119 0.9940 0 0

-1.0000 1.0000 0 0 0 0 1.0000 0

0 0 -1.0000 1.0000 -1.0000 1.0000 0 1.0000

B1 =

-0.5429 0

-0.2525 0

-0.0047 0

20.2483 0

1.9628 0

2.9803 0

0.0200 0

0 0.0326

0 0.0021

0 0.0000

0 -0.0618

0 -0.0045

0 -0.0123

0 -0.0001

0 0

0 0

C1 =

Columns 1 through 13

0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0

Columns 14 through 16

0 1 0

0 0 1

K =

Columns 1 through 8

0.0367 -0.0578 0.0407 0.0634 0.0239 -0.0813 0.1013 0.0485

21.0412 -24.2138 21.6345 4.5472 11.1495 -21.5375 25.0390 22.7474

Columns 9 through 16

-0.0419 0.0319 0.0446 0.0349 -0.0865 0.1026 0.0000 0.0001

-21.8436 20.3135 0.9756 13.2017 -22.4572 25.2658 0.0195 0.0270

L =

0.0925 -0.0000

0.1180 0.0000

0.3752 0.0000

0.0000 0.0568

0.0000 0.1807

0.0000 0.0279

-0.0000 0.2379

0.0694 -0.0000

0.0275 0.0000

0.2971 0.0000

-0.0000 0.0629

-0.0000 0.1964

0.0000 0.0673

0.0000 0.3243

1.6139 0.0000

0.0000 1.6702

5 Моделирование замкнутой системы и оценка качества переходных процессов

Рисунок 5 – Структурная схема системы в виде переменных состояния с учетом запаздывания.


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

AA=[A1 -B1*K;L*C1 A1-B1*K-L*C1];

BB=[B1;zeros(n+l,m)];

CC=[C zeros(l) zeros(l,n+l)];

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

Максимальное возмущение следует принять на уровне 10% номинального значения соответствующих параметров. Допустимое значение урегулированных переменных нужно принять равными 20% номинального значения. Если качество не соответствует нужно сменить весовые матрицы и повторить расчет.

Рисунок 6.Переходные процессы замкнутой системы.


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

6 Преобразование модели регулятора в форму, отвечающую ее реализации в програмном обеспечении

ПИ закон регулирования вычисляется по формулам:

K1=K(:,1:14);

K2=K(:,15:16);

L1=L(1:14,:);

L2=L(15:16,:);

Ar=[Ad-Bd*K1 -Bd*K2-L1 L1; C eye(2)-L2 L2; zeros(2,14) zeros(2) eye(2)];

Br=[zeros(14,2); zeros(2); eye(2)];

Cr=[-K zeros(2)];

Az=[Ad Bd*Cr; Br*C Ar];

Bf=[Bd; zeros(18,2)];

Bz=[zeros(14,2); Br];

Cz=[C zeros(2,18)];

и записать замкнутую систему в вид в котором она будет реализована в программном обеспечении:

x=zeros(14,1);xr=zeros(18,1); u=zeros(2,1);

yy=[]; uu=[];f=[.0010;.0010];z=[0;0];

for i=1:2000,

y=C*x; e=-z+y;

u=Cr*xr; xr=Ar*xr+Br*e;

y=C*x; x=Ad*x+Bd*(u+f);

yy=[yy; y']; uu=[uu; u'];

end

x1=x; xr1=xr; u1=u;

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

 

7 Выбор технических средств реализации системы управления

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

Общая структурная схема рабочей станции изображена на рисунке

Рабочая станция имеет вид:

Рисунок 7 – Схема рабочей станции


Таблица 3

Спецификация технических средств:

Тип К-во Предназначение
1б,2б ADAM 4012 2 Модуль аналогового ввода с датчиков давления, спиртометра, тип входного сигнала: mV, V или mA, диапазон: ±150мВ, ±20мА, ±5В, ±10В
ADAM 4561 1 Преобразователь интерфейса USB в RS-232/422/485
Сапфир-22М-ДА2050 1

Датчик давления в магистрали, верхние пределы измерений: 1.6 МПа, Предел допускаемой основной погрешности, 0,5

0,25; 0,5%, вых. 0-5; 4-20; 5-0; 20-4мА,

Alcolyzer Plus Spirits 1

Спиртомер для крепких спиртных напитков; Диапазон измерения: Спирт: от 35 до 65 об. % (значения отображаются от 0 до 90 об. % спирта, однако, при содержании спирта менее 35 об. % и более 65 об. % точность измерения уменьшается);

Значение рН (опционально): от 0 до 14;

Цвет (опционально): от 0 до 120 EBC;

Плотность (опционально): от 0 до 3 г/см3;

2д,3д ADAM 4069 2 Модуль c релейными выходами, 8 реле с нормально разомкнутым контактом, нагрузочная способность контактов: 250 В/ 5 A для перем. тока, 30 В/ 5 A для пост. тока, время включения 5 мс , время выключения 5,6 мс

2г,2в,

3в,3г

МЭО 40/25-0,25 4 Механизм исполнительный одно-оборотный , номинальный крутя-щий момент 40кгс/м, номинальный ход выходного органа 0,25 оборота за 25с, Напряжение питания 220В. Частота 50Гц
2г,2в КРП-100 2 Клапан регулирующий, Ду=50мм
3г,3в КРП-100 2 Клапан регулирующий, Ду=200мм

Рисунок 8. Функциональная схема автоматизации.


Вывод

В курсовом проекте было выполнено математическое обеспечение АСУТП колонны перегонки коньячного спирта. Были построены характеристики объекта: кривые разгона и при 10-и процентном возмущении. Результаты показали, что качество отвечает требуемому. В результате выбора технического обеспечения: разработана функциональная схема автоматизации, подобрано оборудование для технической реализации данной системы. Разработано программное обеспечение: программа, которая моделирует поведение системы. При тестировании данной программы было показано регулятор работает адекватно.


Литература

1.Стопакевич А.А.Теория систем и системный анализ.Учебник для вузов.-Киев: ВИПОЛ,1996.-200с.

2.Демченко В.А.Автоматизация и моделирование технологичных процессов АЭС и ТЭС.-Одесса:Астроприт,2001.-308с.

3.Стопакевич А.А.Матлаб. Методические указания к лабораторным работам,курсового и дипломного проектирования.-Одесса,2000.-18с.


Приложение 1

 

Текст программы

A=[-1/35 0 0 0 0 0 0 0 0 0 0 0 0 0;

0 -1/129 0 0 0 0 0 0 0 0 0 0 0 0;

0 4/48 -2/48 0 0 0 0 0 0 0 0 0 0 0;

0 0 0 -1/38 0 0 0 0 0 0 0 0 0 0;

0 0 0 4/9 -2/9 0 0 0 0 0 0 0 0 0;

0 0 0 0 0 -1/110 0 0 0 0 0 0 0 0;

0 0 0 0 0 4/134 -2/134 0 0 0 0 0 0 0;

0 0 0 0 0 0 0 -1/13.5 0 0 0 0 0 0;

0 0 0 0 0 0 0 0 -1/98 0 0 0 0 0;

0 0 0 0 0 0 0 0 4/133 -2/133 0 0 0 0;

0 0 0 0 0 0 0 0 0 0 -1/50 0 0 0;

0 0 0 0 0 0 0 0 0 0 4/12 -2/12 0 0;

0 0 0 0 0 0 0 0 0 0 0 0 -1/186 0;

0 0 0 0 0 0 0 0 0 0 0 0 4/150 -2/150];

B=[-42.5/35 0;

-72.5/129 0;

0 0;

1720/38 0;

0 0;

730/110 0;

0 0;

0 0.994/13.5;

0 0.459/98;

0 0;

0 -6.9/50;

0 0;

0 -5.1/186;

0 0];

C=[1 -1 1 0 0 0 0 1 -1 1 0 0 0 0;

0 0 0 -1 1 -1 1 0 0 0 -1 1 -1 1];

D=[0 0;0 0];

dt=0.1/max(abs(eig(A)));

t=0:dt:999;

G=length(t);

[Ad Bd]=c2d(A,B,dt);

y=dstep(Ad,Bd,C,D,1,G);

figure(1)

subplot(2,2,1);plot(y(:,1));grid;ylabel('u1,МПа');title('Razgon u1,1%');

subplot(2,2,3);plot(y(:,2));grid;ylabel('u2,');

%u2

y=dstep(Ad,Bd,C,D,2,G);

subplot(2,2,2);plot(y(:,1));grid;ylabel('y1,М3/с');title('Razgon u2,1%');

subplot(2,2,4);plot(y(:,2));grid;ylabel('y2,кг/с');

A1=[Ad zeros(14,2);C eye(2)];

B1=[Bd; zeros(2)];

C1=[zeros(2,14) eye(2)];

Q2=1e6*[1 0;0 1];

V=C'*Q2*C;

Q=[V zeros(14,2);zeros(2,14) eye(2)];

R=eye(2);

Q1=eye(16);

R1=eye(2);

K=dlqr(A1,B1,Q,R);

L=dlqr(A1',C1',Q1,R1)';

K1=K(:,1:14);

K2=K(:,15:16);

L1=L(1:14,:);

L2=L(15:16,:);

Ar=[Ad-Bd*K1 -Bd*K2-L1 L1; C eye(2)-L2 L2; zeros(2,14) zeros(2) eye(2)];

Br=[zeros(14,2); zeros(2); eye(2)];

Cr=[-K zeros(2)];

Az=[Ad Bd*Cr; Br*C Ar];

Bf=[Bd; zeros(18,2)];

Bz=[zeros(14,2); Br];

Cz=[C zeros(2,18)];

x=zeros(14,1);xr=zeros(18,1); u=zeros(2,1);

yy=[]; uu=[];f=[.0010;.0010];z=[0;0];

for i=1:2000,

y=C*x; e=-z+y;

u=Cr*xr; xr=Ar*xr+Br*e;

y=C*x; x=Ad*x+Bd*(u+f);

yy=[yy; y']; uu=[uu; u'];

end

x1=x; xr1=xr; u1=u;

figure(2)

subplot(2,2,1); plot(yy(:,1));grid;ylabel('y1,MPa');title('Perehod proces Braga ');

subplot(2,2,3); plot(yy(:,2));grid;ylabel('y2,');

subplot(2,2,2); plot(uu(:,1));grid;ylabel('u1,M3/c');title('Perehod proces Par ');

subplot(2,2,4); plot(uu(:,2));grid;ylabel('u2,kg/c');


Еще из раздела Экономико-математическое моделирование:


 Это интересно
 Реклама
 Поиск рефератов
 
 Афоризм
Самый короткий день в году – 1 января: просыпаешься, а за окном уже темнеет…
 Гороскоп
Гороскопы
 Счётчики
bigmir)net TOP 100