Научная статья на тему 'Совместное решение обратной задачи теплопроводности и задачи оптимального проектирования в технологии сварки неплавящимся электродом'

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

CC BY
339
78
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
МОДЕЛИРОВАНИЕ СВАРОЧНЫХ ПРОЦЕССОВ / ПАРАЛЛЕЛЬНЫЕ ВЫЧИСЛЕНИЯ / ЗАДАЧА ОПТИМИЗАЦИИ / ОБРАТНАЯ ЗАДАЧА ТЕПЛОПРОВОДНОСТИ / ЭЛЕКТРОННЫЙ РЕГЛАМЕНТ / ВИРТУАЛЬНОЕ РАБОЧЕЕ МЕСТО / WELDING SIMULATING / PARALLEL COMPUTING / OPTIMIZATION PROBLEM / INVERSE HEAT TRANSFER PROBLEM / DIGITAL REGULATIONS / VIRTUAL WORKPLACE

Аннотация научной статьи по математике, автор научной работы — Кректулева Раиса Алексеевна, Батранин Андрей Викторович

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

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

Похожие темы научных работ по математике , автор научной работы — Кректулева Раиса Алексеевна, Батранин Андрей Викторович

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

The article considers the problems of optimal design in welding technique based on solution of the inverse heat transfer problems applying parallel computing and development of digital regulations for maintaining quality of the obtained process solutions. The authors introduce the solution technique and diagram of parallel search for welding mode optimal parameters. The results obtained by the program «Virtual workplace for welding engineer» developed by the authors are used in the article.

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

УДК 53.09:621.791

СОВМЕСТНОЕ РЕШЕНИЕ ОБРАТНОЙ ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ И ЗАДАЧИ ОПТИМАЛЬНОГО ПРОЕКТИРОВАНИЯ В ТЕХНОЛОГИИ СВАРКИ НЕПЛАВЯЩИМСЯ ЭЛЕКТРОДОМ

Р.А. Кректулева, А.В. Батранин

Томский политехнический университет E-mail: batranin@tpu.ru

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

Ключевые слова:

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

Key words:

Welding simulating, parallel computing, optimization problem, inverse heat transfer problem, digital regulations, virtual workplace.

Введение

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

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

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

В качестве вычислительной базы использовался суперкомпьютерный кластер «СКИФ-политех» [3] с установленным программным обеспечением «Виртуальное рабочее место инженера-сварщика» (ВРМ) [4, 5], разработанным авторами.

Постановка задачи теплопроводности

в технологии сварки

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

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

Математическая постановка задачи моделирования сварочных процессов подробно изложена в работах [6-8]. Задача теплопереноса решается на основе закона сохранения энергии тепловых

процессов и рассматривается пространственное уравнение теплопроводности в виде:

С(Т)Р(Г)^ = div{A(T) • grad(r)} + qv(х y z, t), (1) dt

где С(Т) — удельная теплоемкость, Дж/(кг-К); р(Т) — плотность, кг/м3; Т — температура, К; А(Т) -теплопроводность, Вт/(м-К); qv(x,y,z,t) — суммарная мощность источников тепловой энергии в данной точке пространства, Вт/м3.

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

T (х, y, z,0) = f (х, y, z) (2)

(на практике, как правило, fx,y,z)=const=r0), а также граничными условиями (ГУ), которые в случае дуговых процессов имеют вид:

• вне области действия источника:

дТ

Я — = a(T - To) + ест (Т4 - То4); дп

в области действия источника:

Я — = nUk exp,-kr >),

дп ж

(3)

(4)

где е - степень черноты тела; а=5,669-10-8 Вт/(м2-К4) - постоянная Стефана-Больцмана; п - нормаль к поверхности, м; а - коэффициент теплообмена с окружающей средой, Вт/(м2-К); п - термический КПД; I - ток дуги, А; и - напряжение на дуге, В; к - коэффициент концентрации дуги, 1/м2; г - расстояние от точки до оси дуги, м.

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

x = x0 + Vxt

y = Уо + Vyt,

(5)

V -

где х0, у0 - координаты начальной точки; V скорости движения источника вдоль осей, м/с; / -время сварки, с. Тогда расстояние г из (4) опреде ляется как

Г = 4(x - xo)2 + (У - -Уо)2.

(б)

Следует отметить, что в автоматизированных сварочных процессах движение происходит лишь вдоль одной оси, и это упрощает решаемую задачу (1-6).

Уравнение (1) решается методом конечных разностей [9]. Разностный вид уравнения с учетом нелинейности теплопроводности и неоднородности материала:

C (T, ) ph2

(Т/,,+1 - Tj

+(T,

+(Т+,

j+l,,

- Tk,,, ) Я,

+ (Jk - т

j,,+1/2 + (Т, j^-1 Т j,,

', j+1/2,,

+ (Tik,

) Д, j,,-1/2 +

- Tj ) Я,-1/2,, ■

, - т

(Tj+1 - т j

+(Ti,

j+l,,

) Л+1/2, /,,

)2 Я

- t;* ,, )2 я,,j+1/2,, + (t,

., + 1/2 + (Tikj ,,- 1 - Tik,,,

)Яі-1/2, j,,

)2 Яі, j,,-1/2 +

- Tikj,, )2Яі^,-1/2,, +

+(T+l,j,, - Tk,,, )2 Я.+1/2,,,, + (T-1,j,, - Tikj,,)2 Я,-1/

(7)

где А/ - шаг по времени, с; к - шаг по пространству, м;_Х1,1,1+112—(Х1^+1+Х1^/7, Л^,1-1/2—(Л^,1-1+ и т. д.; Луд1/2—(ЗЛуд1/2)/0Т - скорость изменения теплопроводности, которую находим путем дифференцирования функции Л(Т).

Значения коэффициента теплопроводности Л^ в узлах ячеек соответствуют зависимости Л(Т), определяемой из эксперимента. В общем случае (с учетом структурно-фазовых особенностей материала) экспериментальные данные аппроксимиро-

6

вали полиномом шестой степени Л(Т) = ^ аиГ”.

п=0

Аналогичным образом аппроксимировали экспериментальные данные для теплоемкости

6

С (Т) = ^ЬИГ”. Если в экспериментах структурно-

п=0

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

Р=Ро(1 + рт)-1,

где р0 - плотность материала при нормальных условиях, кг/м3; в - коэффициент теплового расширения, К-1.

Аппроксимацию ГУ при решении пространственных уравнений теплопроводности проводили со вторым порядком точности, что позволяет сохранить передаваемую через границу энергию и обеспечить сходимость решения. С этой целью использован алгоритм аппроксимации ГУ по аналогии с работой [10]. Так, на границе х—0 при заданном тепловом потоке дЧ1 из условия равенства потоков справа и слева от границы, имеем разностное уравнение:

( -

Tlk,,, (2\,,,) + Tk,; До,, - ,

h)

з-

Я1,

Я

-h2 • с(Ток,.,,)p(Tkj.,,)

rpk+1 rpk

то,j,, - То,J,,

At

= о,

(8)

аппроксимирующее ГУ (ГУ 2-го рода) с точностью 0(й2). Аналогично записываются уравнения для остальных границ образца.

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

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

X

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

Рис. 1. Конструктивные элементы сварного шва типа О по ГОСТ 14771-76: е - ширина шва, д ид1 - усиление шва и обратного валика

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

В большинстве случаев для различных материалов и типов сварных соединений, стандарты устанавливают следующие параметры качества: ширина шва е, высота шва или глубина проплавления Н (определяется толщиной изделия и числом проходов), катет шва К (для угловых швов), величина усиления шва g и обратного валика gl. Отметим, что общепризнанной методики выбора параметров режима аргонодуговой сварки пока не предложено, поэтому авторы опираются на методику, разработанную для сварки в защитном газе плавящимся электродом. Согласно ей из приведенного выше набора геометрических параметров вычисляются ширина шва и глубина проплавления (высота шва) [12], которые описывают требуемую изотерму плавления.

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

Согласно формуле (4), в общем случае имеются 4 переменные, которые определяют фактическое значение передаваемой энергии: I, V, п, к. Управляемыми из них являются первые две: ток и напряжение на дуге. Термический КПД и коэффициент концентрации дуги - неуправляемы и определяются экспериментальным путем. Однако в одном и том же технологическом процессе, т. е. при заданном способе сварки и в определенном диапазоне сварочных токов и напряжений, эти величины изменяются в небольших пределах, поэтому в данном случае можно их считать постоянными.

На практике ток и напряжение - величины дискретные и заменяются соответствующими векторами (одномерными массивами) I, Ц, где г-е иу-е значение величин формируются с некоторым шагом в выбранном диапазоне. Имея в виду, что оптимизируемая тепловая энергия 0уд=/(/, V) есть функция двух переменных, переходим к матрице 0у=1-х Ц, которая трансформируется в вектор удельных мощностей источника энергии 0„, где п=г'ху. Таким образом, получается первый управляемый параметр. Другим управляемым параметром процесса служит скорость сварки У, также представляющий собой вектор Ук.

Рассматриваемая задача может быть выражена в терминах задач нелинейного программирования [11]. Исходя из сказанного, получаем двухпараметрическую нелинейную_ задачу оптимизации необходимой энергии Д0п,Ук)^шт со следую -щими ограничениями:

О < О < О ;

£-т1п ^~-п 2^-тах ’ 0 < V, < Ктах; т Л) = е;

/2(Оп Л) = А. (9)

Здесь Е(0п,¥к) - целевая функция; 0п и Ук -искомые векторы. 0шп, йшж, Уши будут определяться возможностями сварочного оборудования и условиями физического существования процесса. Функции /1 и /2 - нелинейны и определяются из результатов компьютерного эксперимента. Отметим, что получение значений аргументов функций /1 и /2, отвечающих условию (9), представляет собой решение обратной задачи. Таким образом, решение задачи оптими_зац_ии есть поиск минимума целевой функции Е( 0п, Ук) на матрице размером пхк с учетом системы ограничений (9).

Методология использования «Виртуального

рабочего места» для решения поставленных задач

«Виртуальное рабочее место» (ВРМ) разрабатывалось с использованием представленной математической модели (1-6) для решения, прежде всего, прямых задач технологии сварки, а именно, для получения параметров сварных соединений по заданным технологическим режимам.

Функции /1 и /2 из (9) не могут быть обращены в большинстве случаев в связи со сложностью по-ставленн_ой за_дачи, поэтому для получения аргументов 0п и Ук необходимо решать набор прямых зада_ последовательно перебирая элементы векторов 0п и Ук. Чисто технически, данный подход затратен по времени, если весь набор решений получать последовательно. Поэтому представляется целесообразным параллельное решение прямых задач с использованием вычислительного кластера.

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

ный кластер «СКИФ-политех», часть вычислительных узлов которого работает под управлением Windows Server 2008 HPC, позволяет работать с несколькими копиями ВРМ, функционирующими в системе Matlab.

Рассмотрим пример решения задачи оптимального проектирования технологического процесса сварки с применением ВРМ. Допустим, необходимо спроектировать сварное соединение типа С7 по ГОСТ 14771-76 (рис. 1) из низкоуглеродистой стали для деталей толщиной 6 мм. Моделируется сварка в защитном газе неплавящимся электродом без присадочной проволоки. Соединение без зазора, сварка выполняется в два прохода с лицевой и обратной стороны. Так как швы симметричны и выполняются на одном и том же режиме, то рассмотрим лишь процесс выполнения одного шва. Согласно ГОСТ, имеем ширину шва не более 10 мм. Чтобы гарантировать проплавление, зададим высоту шва для одного прохода не менее 3 мм.

Численный эксперимент проводился на образце размером 80x50x6 мм с прямолинейно движущимся источником по центру образца (рис. 2).

дач образуют диагональ матрицы пхк—и дальнейший поиск минимума функции Е(0п,¥к) производится на этой диагонали. В таблице представлены некоторые варианты технологических процессов, которые были исследованы на поиск минимума функции Е(0п,¥к). Во всех случаях в формуле (4) принято: к=7 см-2, КПД 90 %, граничные условия - адиабатические (теплообмен отсутствует).

Таблица. Параметры процесса сварки (численный эксперимент)

Параметр Режим № 1 Режим № 2 Режим № 3 Режим № 4

Ток, А 180 220 310 350

Напряжение, В 11 13 17 18

Скорость сварки, мм/с 2,0 3,6 7,8 9,4

Мощность источника нагрева, кВт 1,78 2,57 4,74 5,67

Время сварки образца, с 40,0 22,2 10,3 8,5

Затраченная энергия, кДж 71,3 57,2 48,7 48,3

Погонная энергия, 106 Дж/м 0,89 0,72 0,61 0,60

Рис. 2. Модель экспериментального образца

Для проведения эксперимента устанавливались контрольные точки, которые расположены в характерных областях: на поверхности образца на расстоянии 5 мм от оси траектории движения источника нагрева, чтобы отслеживать ширину шва, и на глубине 4 мм от поверхности на оси траектории, чтобы контролировать глубину проплавления. Схема расположения контрольных точек в поперечном относительно движения источника нагрева сечении приведена на рис. 3.

Рис. 3. Расположение контрольных точек в поперечном сечении образца

Получение решений обратных задач требует пхк решений прямых задач. Решения обратных за-

Термограммы в контрольных точках для всех режимов сварки приведены на рис. 4. Температура плавления низкоуглеродистой стали 1800 К.

Из рис. 4 следует, что температуры в контрольных точках примерно равны температуре плавления: условие (9) выполняется. Исключение составляет контрольная точка на поверхности при сварке на режиме № 4. Здесь температура выше на 200 К, что говорит о превышении допустимой ширины шва: условие (9) не выполняется. Допустимым режимом сварки можно считать режим № 3, который не приводит к заметному перегреву и превышению ширины шва. Рис. 5 иллюстрирует решение задачи оптимального проектирования.

Таким образом, оптимальным будет режим № 3. Ему соответствует минимальное время сварки и минимальная затраченная энергия при одновременном выполнении условий (9).

Принципиальная схема параллельного решения задачи приведена на рис. 6.

Исходные данные: материал и геометрия образца, пределы изменений параметров режима, требуемые свойства сварного шва и пр. - поступают в программу-диспетчер, которая распределяет подзадачи между узлами кластера - решателями. Подзадачи представляют собой расчет отдельных процессов, параметры режима которых выбраны из диапазона матрицы пхк. По окончании расчета подзадачи решатель передает полученные данные диспетчеру. Конечным этапом - получением решения оптимизационной -адачи - является поиск минимума функции Е (0п,Ук), как показано выше. Результаты решения вместе с исходными данными заносятся в базу данных для дальнейшего использования, например, для пополнения электронного регламента на обеспечение качества технологических решений.

Режим № 1

Режим № 3

Режим № 2

Режим № 4

Рис. 4. Термические циклы в глубине образца (-----) и на границе шва (-----)

Е, кДж 80^.

Рис. 5. Решение задачи оптимизации на основе численных экспериментов

Исходные данные

Диспетчер - * Решение

Система очередей База данных

г

л л

5 §

со Э со

0 0)

0. О.

Рис. 6. Схема параллельного поиска оптимального решения

Электронный регламент

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

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

• внести критерии качества продукции;

• отразить способ оптимизации и управления технологическими процессами;

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

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

В целом разрабатываемый электронный регламент должен соответствовать ГОСТ Р 52294-2004 «Информационная технология. Управление организацией. Электронный регламент административной и служебной деятельности. Основные положения».

СПИСОК ЛИТЕРАТУРЫ

1. Валдайцева Е.А., Туричин ГА., Норман Е.А., Земляков Е.В., Малкин П.Е. Методика выбора режимов лазерной сварки средствами LaseгCAD // Лучевые технологии и применение лазеров: Труды V Междунар. научно-техн. конф. - г. Санкт-Петербург, 23-28 сентября 2006. - СПб.: Изд-во СПбГПУ, 2006. - С. 307-316.

2. Ерофеев В.А. Решение задач оптимизации технологии на основе компьютерного моделирования процесса сварки // Сварочное производство. - 2003. - № 7. - С. 19-26.

3. Программирование, компиляция и запуск программ с использованием вычислительного кластера «СКИФ-политех». Материалы Лаборатории мультифизического моделирования на базе вычислительного кластера «СКИФ-политех» ТПУ // Томский политехнический университет. 2011. ТОЬ: Шр://с1и-ster.tpu.ru (дата обращения: 10.02.2011).

4. Батранин А.В., Кректулева Р.А., Черепанов Р.О. Виртуальное рабочее место инженера-сварщика (ВРМ). Свидетельство РФ о государственной регистрации программы для ЭВМ № 2010616487 (дата регистрации: 30.09.2010).

5. Кректулева Р. А., Батранин А. В. Разработка виртуального рабочего места для подготовки инженеров-сварщиков // Современные проблемы машиностроения: Матер. IV Междунар. на-учно-техн. конф. - г. Томск, 26-28 ноября 2008. - Томск: Изд-во ТПУ, 2008. - С. 319-323.

Выводы

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

2. Поставлена и решена задача оптимизации сварочных процессов, основанная на решении обратных задач теплопроводности.

3. Приведена схема параллельного решения задачи с использованием вычислительного кластера и разработанного программного обеспечения. Рекомендовано к публикации Оргкомитетом V Междунар. научн-техн. конф. «Современные проблемы машиностроения», Томск, ТПУ, 23—26 ноября 2010 г.

6. Кректулева РА., Никифоров Н.И., Сухинин Г.К., Бежин О.Н., Губенко Л.В. Результаты компьютерных и натурных экспериментов по высокоскоростной кислородной резке металла // Aвтоматическая сварка. - 2000. - N° 5. - С. 21-24.

7. Кректулева РА., Сараев Ю.Н., Косяков ВА Математическое моделирование технологических процессов импульсной аргонодуговой сварки неплавящимся электродом // Сварочное производство. - 1997. - М 4. - С. 2-4.

8. Кректулева РА., Батранин A^., Бежин О.Н. Применение программного обеспечения Meza для оценки дефектности сварных соединений на стадии проектирования // Сварка и диагностика. - 2009. - М 2. - C. З6-42.

9. Самарский A.A., Вабищевич П.Н. Вычислительная теплопередача. - М.: Едиториал УРСС, 200З. - 784 с.

10. Дульнев Г.Н., Парфенов В.Г., Сигалов A^. Применение ЭВМ для решения задач теплообмена. - М.: Высшая школа, 1990. -207 с.

11. Карманов В.Г. Математическое программирование. - М.: Физ-матлит, 2008. - 264 с.

12. Aкулов A.K, Бельчук ГА., Демянцевич В.П. Технология и оборудование сварки плавлением. - М.: Машиностроение, 1977. -4З2 с.

Поступила 10.02.2011 г.

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