Идентификация распределенных параметров внешнего теплообмена для нелинейных граничных условий
Збережено в:
| Опубліковано в: : | Труды Института прикладной математики и механики НАН Украины |
|---|---|
| Дата: | 2008 |
| Автор: | |
| Формат: | Стаття |
| Мова: | Russian |
| Опубліковано: |
Інститут прикладної математики і механіки НАН України
2008
|
| Онлайн доступ: | https://nasplib.isofts.kiev.ua/handle/123456789/20015 |
| Теги: |
Додати тег
Немає тегів, Будьте першим, хто поставить тег для цього запису!
|
| Назва журналу: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Цитувати: | Идентификация распределенных параметров внешнего теплообмена для нелинейных граничных условий / О.В. Литовченко // Труды Института прикладной математики и механики НАН Украины. — Донецьк: ІПММ НАН України, 2008. — Т. 17. — С. 109-118. — Бібліогр.: 6 назв. — рос. |
Репозитарії
Digital Library of Periodicals of National Academy of Sciences of Ukraine| id |
nasplib_isofts_kiev_ua-123456789-20015 |
|---|---|
| record_format |
dspace |
| spelling |
Литовченко, О.В. 2011-05-20T07:41:52Z 2011-05-20T07:41:52Z 2008 Идентификация распределенных параметров внешнего теплообмена для нелинейных граничных условий / О.В. Литовченко // Труды Института прикладной математики и механики НАН Украины. — Донецьк: ІПММ НАН України, 2008. — Т. 17. — С. 109-118. — Бібліогр.: 6 назв. — рос. 1683-4720 https://nasplib.isofts.kiev.ua/handle/123456789/20015 51-74:536.2 ru Інститут прикладної математики і механіки НАН України Труды Института прикладной математики и механики НАН Украины Идентификация распределенных параметров внешнего теплообмена для нелинейных граничных условий Article published earlier |
| institution |
Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| collection |
DSpace DC |
| title |
Идентификация распределенных параметров внешнего теплообмена для нелинейных граничных условий |
| spellingShingle |
Идентификация распределенных параметров внешнего теплообмена для нелинейных граничных условий Литовченко, О.В. |
| title_short |
Идентификация распределенных параметров внешнего теплообмена для нелинейных граничных условий |
| title_full |
Идентификация распределенных параметров внешнего теплообмена для нелинейных граничных условий |
| title_fullStr |
Идентификация распределенных параметров внешнего теплообмена для нелинейных граничных условий |
| title_full_unstemmed |
Идентификация распределенных параметров внешнего теплообмена для нелинейных граничных условий |
| title_sort |
идентификация распределенных параметров внешнего теплообмена для нелинейных граничных условий |
| author |
Литовченко, О.В. |
| author_facet |
Литовченко, О.В. |
| publishDate |
2008 |
| language |
Russian |
| container_title |
Труды Института прикладной математики и механики НАН Украины |
| publisher |
Інститут прикладної математики і механіки НАН України |
| format |
Article |
| issn |
1683-4720 |
| url |
https://nasplib.isofts.kiev.ua/handle/123456789/20015 |
| citation_txt |
Идентификация распределенных параметров внешнего теплообмена для нелинейных граничных условий / О.В. Литовченко // Труды Института прикладной математики и механики НАН Украины. — Донецьк: ІПММ НАН України, 2008. — Т. 17. — С. 109-118. — Бібліогр.: 6 назв. — рос. |
| work_keys_str_mv |
AT litovčenkoov identifikaciâraspredelennyhparametrovvnešnegoteploobmenadlânelineinyhgraničnyhuslovii |
| first_indexed |
2025-11-27T07:52:35Z |
| last_indexed |
2025-11-27T07:52:35Z |
| _version_ |
1850804281928581120 |
| fulltext |
ISSN 1683-4720 Труды ИПММ НАН Украины. 2008. Том 17
УДК 51-74:536.2
c©2008. О.В. Литовченко
ИДЕНТИФИКАЦИЯ РАСПРЕДЕЛЕННЫХ ПАРАМЕТРОВ
ВНЕШНЕГО ТЕПЛООБМЕНА ДЛЯ НЕЛИНЕЙНЫХ
ГРАНИЧНЫХ УСЛОВИЙ
Исследуется возможность идентификации распределенных параметров модели теплофизического
процесса с помощью метода наименьших квадратов. Рассмотрены случаи идентификации функции
параметра лучистого теплообмена и функции параметра конвективного теплообмена, имеющей два
локальных максимума. Представлен алгоритм идентификации, приведены и проанализированы
данные численных расчётов по нахождению распределенных по времени параметров.
1. Введение. Задача идентификации распределенного параметра относится к
обратным задачам и является некорректно поставленной в классическом смысле.
К настоящему времени наиболее популярным методом решения такого класса за-
дач является метод регуляризации Тихонова. Недостатком метода регуляризации
является большой объем вычислений, связанный с процедурой поиска нужного зна-
чения параметра регуляризации, а также с овражностью регуляризирующего функ-
ционала [4]. Другой существенный недостаток метода связан с самой идеей регуля-
ризации: сглаживания решения в пределах погрешности измерений. Чем больше
погрешность, тем можно получить более гладкую кривую, но при этом возрастает
опасность получения хотя и более плавной кривой, но все в большей степени откло-
няющейся от истинной [1, 3]. В этом смысле метод регуляризации принципиально
отличается от ряда методов, которые обеспечивают получение достаточно точного
решения при росте погрешности в исходных данных. К таким методам относится,
например, метод наименьших квадратов.
В представленной работе для идентификации распределенных параметров мо-
дели теплофизического процесса предлагается применить метод наименьших квад-
ратов.
2. Идентификация параметра лучистого теплообмена. Сформулируем об-
ратную задачу нахождения коэффициента лучистого теплообмена для одномерно-
го уравнения теплопроводности с нелинейными граничными условиями 3-го рода.
Ставится задача нахождения σ = σ(τ) как функции, зависящей от времени в виде
аппроксимирующего полинома.
Математическая модель.
Математическая модель процесса выглядит следующим образом
∂ T
∂τ
= a
∂2 T (τ, x)
∂x2
, 0 ≤ x ≤ ` (1)
с нелинейными граничными условиями третьего рода
−λ
∂T (τ, x)
∂x
∣∣∣∣
x=0
= σ1
[
T 4
gr.(τ)− T 4(τ, o)
]
, (2)
109
О.В. Литовченко
λ
∂T (τ, x)
∂x
∣∣∣∣
x=l
= σ2
[
T 4
gr.(τ)− T 4(τ, l)
]
, (3)
и начальным условием
T (0, x) = t0(x), (4)
где T (τ, x) – температура тела, Tgr.(τ) – температура греющей среды, λ – коэффици-
ент теплопроводности среды, σ1, σ2 – коэффициенты лучистой теплоотдачи сверху
и снизу.
Для решения обратных задач обычно требуются знания о температуре нагрева-
емого тела [5].
Для простоты изложения метода предполагаем нагрев симметричным, поэтому
σ1 = σ2 = σ. В этом случае минимальный объем необходимой для решения задачи
информации соответствует измерению температуры тела в какой-либо одной точке
одномерной области o ≤ x ≤ l. Предположим, известна температура тела на границе
с внешней средой:
T (τ, 0) = f(τ). (5)
Задача состоит в нахождении σ = σ(τ) как функции, зависящей от времени. Пред-
полагая функцию σ(τ) непрерывной, с целью аппроксимации искомой функции вос-
пользуемся полиномом n-ой степени
σ(τ) = a0 + a1 · τ + ... + an · τn, (6)
степень полинома будет определяться по принципу невязки.
Решение обратной задачи.
Решив задачу Дирихле, получим температуры T (τi, xj), которые будут необхо-
димы для вычисления производной в граничном условии. Для решения задачи Ди-
рихле используем метод конечных разностей и метод прямой итерации.
В области 0 ≤ x ≤ l, 0 ≤ τ ≤ τ̃ введем равномерную сетку ωi,j , т.е. будем рас-
сматривать не T (τ, x), а T (i∆τ, j∆x), где ∆x = l
n , i и j шаги по времени и по
координате. Условие сходимости метода 0 ≤ a∆τ
∆x ≤ 1
2 .
Применив явную конечно-разностную схему, получим представление уравнения теп-
лопроводности в виде:
Ti+1,j = c1Ti,j−1 + (1− 2 · c1)Ti,j + c1Ti,j+1,
где c1 = ∆τa
∆x2 .
Теперь в граничном условии λ∂T
∂x
∣∣
x=0
= σ(τ)
[
T 4
gr.(τ)− T 4(τ, 0)
]
неизвестным яв-
ляется только величина лучистого теплообмена σ(τ).
Имеем систему уравнений граничных условий в r моментах времени:
λ
∆x
(T (τi, x1)− T (τi, x0)) = σi[T 4
gr. − T 4(τi, x0)], i = 1, r, (7)
где σk – это значение полинома в момент времени k и σk = a0+a1·k·∆τ+...+ankn∆τn,
k = 1, r, r >> n;
110
Идентификация распределенных параметров
или в матричной форме
λ
∆x
T (τ1, x1)− T (τ1, x0)
T (τ2, x1)− T (τ2, x0)
...
...
T (τr, x1)− T (τr, x0)
=
σ1 · (T 4
gr. − T 4(τ1, x0))
σ2 · (T 4
gr. − T 4(τ2, x0))
...
...
σr · (T 4
gr. − T 4(τr, x0))
=
=
σ1 0 ... 0
0 σ2 ... 0
...
0 0 ... σr
·
T 4
gr. − T 4(τ1, x0)
T 4
gr. − T 4(τ2, x0)
...
...
T 4
gr. − T 4(τr, x0)
. (8)
Выделим компоненты неизвестных параметров полинома, представленного в мат-
ричной форме
σ1 0 ... 0
0 σ2 ... 0
...
0 0 ... σr
=
a0 0 ... 0
0 a0 ... 0
...
0 0 ... a0
+
+
a1 0 ... 0
0 a1 ... 0
...
0 0 ... a1
·
1∆τ 0 ... 0
0 2∆τ ... 0
...
0 0 ... r∆τ
+ ...+
+
an 0 ... 0
0 an ... 0
...
0 0 ... an
·
1∆τn 0 ... 0
0 2n∆τn ... 0
...
0 0 ... rn∆τn
=
= a0 ·
1 0 ... 0
0 1 ... 0
...
0 0 ... 1
+ a1 ·
1∆τ 0 ... 0
0 2∆τ ... 0
...
0 0 ... r∆τ
+ ...+
+an ·
1∆τn 0 ... 0
0 2n∆τn ... 0
...
0 0 ... rn∆τn
.
Для удобства записи преобразований введем такие обозначения
A0 =
1 0 ... 0
0 1 ... 0
...
0 0 ... 1
, A1 =
1∆τ 0 ... 0
0 2∆τ ... 0
...
0 0 ... r∆τ
,
111
О.В. Литовченко
An =
1∆τn 0 ... 0
0 2n∆τn ... 0
...
0 0 ... rn∆τn
, P = λ
∆x
T (τ1, x1)− T (τ1, x0)
T (τ2, x1)− T (τ2, x0)
...
...
T (τr, x1)− T (τr, x0)
,
σ̃ =
σ1 0 ... 0
0 σ2 ... 0.
...
0 0 ... σr
, H =
T 4
gr − T 4(τ1, x0)
T 4
gr − T 4(τ2, x0)
...
...
T 4
gr − T 4(τr, x0)
.
Pk и Hk−k-тые элементы столбцов, Ai
kk−k-й элемент главной диагонали матрицы
Ai.
Тогда уравнение (8) в матричной форме будет иметь вид
P = σ̃H. (9)
Метод наименьших квадратов.
Для нахождения функции методом МНК введем меру отклонения в виде суммы
квадратов разности измеренных температур от расчетных по модели (1)–(6).
S(σ) =
r∑
k=1
(Pk − σkHk)2. (10)
Используя введенные обозначения параметр σ можно записать как
σ(τ) = a0 ·A0 + +a1 ·A1 + ... + an ·An,
тогда квадратичный функционал будет иметь вид
S(σ) =
r∑
k=1
(Pk − σkHk)2 = (P − (a0 ·A0 + a1 ·A1 + ... + an ·An) ·H)T ·
·(P − (a0 ·A0 + a1 ·A1 + ... + anAn) ·H) = P T P − 2
n∑
k=0
akH
T AkP+
+
n∑
i=0
n∑
j=0
aiajH
T Ai+jH.
Задача сводится к нахождению параметров полинома, которые бы минимизиро-
вали функционал S(σ). Пусть, например, минимум достигается при σ′ = σ(τ ′). Тогда
в этой точке должно быть выполнено необходимое условие минимума функционала
S
dS(σ′)
dam
= 0, m = 0, n; (11)
Вычислим векторы первых производных суммы квадратов отклонений по компо-
нентам a0, a1, ..., an :
112
Идентификация распределенных параметров
∂S
∂aj
=
∂[(P−σ̃H)T (P−σ̃H)]
∂aj
= 2a0H
T AjH + 2a1H
T Aj+1H + ...+
+2anHT Aj+nH − 2HT AjP = 0, где j = 0, n.
Получили систему из n+1 линейных алгебраических уравнений с n+1 неизвестными
n∑
m=0
amHT AmH = HT A0P,
n∑
m=0
amHT Am+1H = HT A1P,
. . .
n∑
m=0
amHT Am+nH = HT AnP.
(12)
Из структуры строк матрицы этой системы видна их линейная независимость,
следовательно, детерминант не обращается в ноль и решение системы единственное.
Из системы (12) получаем вектор неизвестных компонент параметра σ(τ)
a0
a1
.
.
.
an
=
HT A0H HT A1H . . . HT AnH
HT A1H HT A2H . . . HT An+1H
. . . . . . . . . . . .
HT AnH HT An+1H . . . HT A2nH
−1
·
HT A0P
HT A1P
. . .
HT AnP
. (13)
Вычислительный эксперимент.
Для проверки работоспособности полученных соотношений была выполнена про-
граммная реализация алгоритма на языке C и проведен численный расчет по нахож-
дению функции σ(τ). В качестве тестовой функции, т.е. истинного распределения
значений параметра, взята функция
σt(τ) · 108 = 0.6 · x4 − 3.02 · x3 + 3.65 · x2 + 0.57 · x + 2.
На известные значения температуры поверхности тела был нанесен шум ξ, ими-
тирующий погрешность измерений [6], распределенный по нормальному закону с
математическим ожиданием равным нулю.
Тогда из соотношений (13) и (6) получили аппроксимирующие полиномы с пер-
вой по четвертую степень (σi, i = 0, 4) (см. рис.1)
Чтобы убедиться в качестве полученных оценок, найдем среднеквадратическое
отклонение (СКО), характеризующее погрешность для накладываемого шума ξ, от-
клонения аппроксимирующего полинома σ(τ) от тестовой функции и для отклоне-
ния расчетной температуры на границе тела от исходных данных.
Расчеты погрешностей СКО найденных решений и отклонения измеренных тем-
ператур от расчетных при погрешности измерения СКОξ = 5.960C приведены в
таблице 1.
113
О.В. Литовченко
Рис. 1. Графики полученных аппроксимирующих полиномов для тестовой функции σt(τ)
Таблица 1. Погрешности найденных решений.
Степень аппроксими-
рующего полинома
σ(τ)
0 1 2 3 4
СКОσ, Вт/м2 ·K4 0.63 0.43 0.14 0.08 0.07
СКОT , 0C 46.40 23.71 21.42 21.01 21.01
Результаты вычислений показывают, что приемлемая точность полученного ре-
шения достигается начиная с полинома 3-ей степени.
3. Идентификация параметра конвективного теплообмена. Рассмотрим
задачу нахождения распределенного во времени параметра конвективного теплооб-
мена α(τ) в граничном условии третьего рода
−λ
∂T (τ, x)
∂x
∣∣∣∣
x=0
= α1 [Tgr.(τ)− T (τ, o)] , (14)
λ
∂T (τ, x)
∂x
∣∣∣∣
x=l
= α2 [Tgr.(τ)− T (τ, l)] , (15)
модели теплового процесса (1), (4), (5), (14), (15). Для простоты предположим α1 =
α2 = α. Для аппроксимации искомой функции воспользуемся полиномом n-ой сте-
пени
α(τ) = a0 + a1 · τ + ... + an · τn, (16)
степень полинома будет определяться по принципу невязки.
Система (7) для случая аппроксимации параметра конвективного теплообмена
принимает следующий вид
λ
∆x
(T (τi, x1)− T (τi, x0)) = αi[Tgr. − T (τi, x0)], i = 1, r, (17)
114
Идентификация распределенных параметров
где αk – это значение полинома в момент времени k и αk = a0+a1·k·∆τ+...+ankn∆τn,
k = 1, r, r >> n.
Уравнение (9) записывается следующим образом
P = α̃H, (18)
где α̃ =
α1 0 ... 0
0 α2 ... 0.
...
0 0 ... αr
, H =
Tgr − T (τ1, x0)
Tgr − T (τ2, x0)
...
...
Tgr − T (τr, x0)
.
Разрешив систему (13), можем найти неизвестные коэффициенты аппроксими-
рующего полинома (16).
Рассмотрим случай идентификации функции α(τ), которая имеет 2 точки ло-
кального максимума на рассматриваемом промежутке времени. Такой вид параметр
конвективного теплообмена может принимать, когда процесс нагрева проводится в
двух рабочих зонах. В качестве тестовой функции была взята
αt(τ) = −224.94 · x4 + 914.52 · x3 − 1165.2 · x2 + 501.81 · x + 30.
После зашумления тестовых значений температуры на границе тела и использова-
ния предложенного метода идентификации из соотношений (13) и (6) для случая
аппроксимации константой, функциями первой-четвертой степеней были получены
следующие результаты (см. Рис.2).
Рис. 2. Графики полученных аппроксимирующих полиномов для тестовой функции αt(τ)
Численные оценки погрешностей найденных решений и отклонения измеренных
температур от расчетных представлены в таблице 2.
115
О.В. Литовченко
Таблица 2. Погрешности найденных решений
Степень аппроксими-
рующего полинома
α(τ)
0 1 2 3 4
СКОα, Вт/м2 ·0 C 26.37 22.52 24.54 27.64 2.46
СКОT ,0 C 29.31 29.04 29.43 24.74 20.55
Полиномы нулевой, первой и второй степеней фактически сгладили экстремумы
тестовой зависимости. Полином 3-ей степени не может в полной мере воспроизвести
искомую функцию. А в случае, когда полином на качественном уровне совпадает с
тестовой функцией, то МНК дает достаточно высокую точность решения задачи.
4. Достаточное условие минимума. Для того чтобы убедиться, что уравне-
ние (13) определяет минимум, достаточно убедиться, что матрица вторых производ-
ных (матрица Гессе) полуположительно определена [2].
G =
HT A0H HT A1H . . . HT AnH
HT A1H HT A2H . . . HT An+1H
. . . . . . . . . . . .
HT AnH HT An+1H . . . HT A2nH
. (19)
Для главных миноров этой матрицы введем такие обозначения: G1 – минор размер-
ностью 1× 1, G2 – минор размерностью 2× 2, Gi – минор размерностью i× i.
Рассмотрим матрицу G: эта матрица не зависит от a0, a1, ..., an, симметрична
и все ее элементы положительны, каждый из них является суммой r слагаемых.
Например,
HT A0H =
r∑
k=1
H2
kA0
kk, HT A1H =
r∑
k=1
H2
kA1
kk, HT AiH =
r∑
k=1
H2
kAi
kk.
Заметим, что матрица представима в виде произведения двух матриц, размер-
ностями n + 1× r и r × n + 1. Причем эти матрицы будут транспонированными по
отношению друг к другу.
Т.е. матрицу (19) можно записать в таком виде
G =
H1A
0
11 H2A
0
22 . . . HrA
0
rr
H1A
1
11 H2A
1
22 . . . HrA
1
rr
. . . . . . . . . . . .
H1A
n
11 H2A
n
22 . . . HrA
n
rr
·
H1A
0
11 H1A
0
11 . . . H1A
0
11
H2A
1
22 H2A
1
22 . . . H2A
1
22
. . . . . . . . . . . .
HrA
n
rr HrA
n
rr . . . HrA
n
rr
. (20)
Обозначим матрицы-множители как Q(n+1× r) и соответственно QT (r×n+1).
Заметим, что любой из главных миноров Gi (i = 1, n + 1) является произведе-
нием первых i строк матрицы Q и первых i столбцов матрицы QT Gi = Qi ·QT
i .
А детерминант Gi равен сумме квадратов детерминантов всех i× i матриц, которые
можно составить из столбцов Qi. Т.е. detGi =
ci
r∑
k=1
(detDk)2, где Dk (k = 1, ci
r) –
описанные i× i матрицы i = 1, n + 1.
116
Идентификация распределенных параметров
Из сказанного выше следует, что detGi > 0 ∀i = 1, n + 1.
По критерию Сильвестра для положительно определенных матриц, из того, что
все главные миноры матрицы Гессе положительны, получаем положительную опре-
деленность матрицы G.
Т.о. мы показали, что условия (11) определяют глобальный минимум функцио-
нала S. И из полученного соотношения (13) можем найти вектор параметра аппрок-
симирующего полинома (6).
5. Свойства оценки параметра. В обратной задаче теплопроводности наря-
ду с температурой поверхности тела используется еще ряд измеряемых величин.
Это такие величины, как: время, координата датчика, толщина образца. Предпо-
лагается, что все они, за исключением температуры, известны точно. А измерения
температуры являются основным источником погрешностей и неопределенностей, и
ее погрешность предполагается белым шумом.
Оценки параметров, получаемые по методу МНК, при условии выполнения пред-
посылок относительно случайных погрешностей наблюдений, будут обладать следу-
ющими свойствами:
1. оценки параметров являются несмещенными, т.е. математическое ожидание
оценок параметров равно истинному значению параметров. Данное свойство яв-
ляется логическим следствием второго предположения о характере погрешности.
Несмещенность означает, что выборочные оценки параметров концентрируются во-
круг неизвестных истинных параметров;
2. оценки состоятельны, иначе говоря, дисперсия оценки параметра стремится к
нулю с возрастанием числа наблюдений r;
3. оценки являются эффективными в том смысле, что они имеют минимальную
дисперсию по сравнению с любыми другими оценками этого параметра.
Если предположение 3 или 4 нарушено, то свойство несмещенности и состоятель-
ности оценок сохраняется, однако оценки оказываются менее эффективными, чем в
случае, когда эти допущения соблюдаются.
Совершенно очевидно, что не безразлично, какими свойствами обладает оценка.
Что касается свойства несмещенности, то оно является необходимым. В самом деле,
смещенные оценки априори дают неверное положение кривой в пространстве неза-
висимых переменных. Свойство состоятельности означает, что при увеличении объ-
ема наблюдения оценки параметров становятся более надежными в вероятностном
смысле, т.е. с ростом r оценки все плотнее концентрируются вокруг истинных неиз-
вестных значений параметров. Свойство эффективности, в общем, является наибо-
лее важным, поскольку оно определяет степень возможной ошибки прогноза.
Основные результаты и выводы. В статье предложен эффективный метод
идентификации распределенных параметров внешнего теплообмена в линейных и
нелинейных граничных условиях задачи теплопроводности. Искомое решение ап-
проксимируется полиномом n-ой степени, коэффициенты которого отыскиваются
методом наименьших квадратов. Разработан алгоритм идентификации, который
программно реализован на языке С. Численные исследования наглядно показали,
что предложенный метод успешно решает поставленную задачу и обладает рядом
117
О.В. Литовченко
преимуществ по сравнению с известными методами решения обратных задач на ос-
нове идеи регуляризации.
1. Алифанов О.М. О методах решения некорректных обратных задач // Инженерно-физический
журнал. – 1983, т.45. – C.742-752.
2. Гантмахер Ф.Р. Теория матриц. М.: Наука, 1966, 576с.
3. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М.: Наука, 1974, C.224
4. Ткаченко В.Н. Моделирование тепловых процессов в автоматизированных системах обработки
информации // Вiсник Донецького нацiонального унiверситету, Серiя А. Природничi науки. –
2002, №2, C.379-383.
5. Мацевитый Ю.М. Обратные задачи теплопроводности. В 2-х т. : Т.2. Приложения. – НАН
Украины, Институт проблем машиностроения. – Киев: Наукова думка, 2003.
6. Коздоба Л.А., Круковский П.Г. Методы решения обратных задач теплопереноса. – Киев: Нау-
кова думка. – 1982. – 385c.
Ин-т прикл. математики и механики НАН Украины, Донецк
lit.ov@i.ua
Получено 04.12.08
118
содержание
Том 17
Донецк, 2008
Основан в 1997г.
|