Научная статья на тему 'Двухуровневая модель статического деформирования конструкций из армированных материалов'

Двухуровневая модель статического деформирования конструкций из армированных материалов Текст научной статьи по специальности «Математика»

CC BY
81
11
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
ИДЕНТИФИКАЦИЯ / IDENTIFICATION / МАТЕМАТИЧЕСКАЯ МОДЕЛЬ / MATHEMATICAL MODEL / КОМПОЗИТНЫЕ МАТЕРИАЛЫ / COMPOSITE MATERIALS / АРМИРОВАННЫЕ МАТЕРИАЛЫ / REINFORCED MATERIALS / НЕРАЗРУШАЮЩИЙ КОНТРОЛЬ / NON-DESTRUCTIVE TESTING

Аннотация научной статьи по математике, автор научной работы — Казначеева Ольга Константиновна, Полинко Юлия Викторовна

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

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

Похожие темы научных работ по математике , автор научной работы — Казначеева Ольга Константиновна, Полинко Юлия Викторовна

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

A two-level model of static deformation of constructions made of reinforced materials

A two-level model consisting of the final element model of static deformation (the basic model) and the reduced model with explicit expressions for response functions is offered. The model represents an inhomogeneous design, with its parameters of elasticity and rigidity representing Piecewise constant functions of coordinates, with external actions being given functions of coordinates and the response being the fields of transitions, strains and stresses. The mathematical model looks like a variation boundary-value problem with variable parametres. For practical calculations it is substituted by the discrete analogue i. e. by the final-element model. The increase of calculation profitability is obtained by introducing the reduced model in the form of explicit approximation of response function dependence on structural parametres. The reduction algorithm of the discrete model is based on the factoral computing experiment. The obtaied reduced model is used for minimization of the quality estimation criterion and for deriving interval estimations of parametres.

Текст научной работы на тему «Двухуровневая модель статического деформирования конструкций из армированных материалов»

УДК 624.04 DOI: 10.17213/0321-2653-2015-4-112-116

ДВУХУРОВНЕВАЯ МОДЕЛЬ СТАТИЧЕСКОГО ДЕФОРМИРОВАНИЯ КОНСТРУКЦИЙ ИЗ АРМИРОВАННЫХ МАТЕРИАЛОВ

A TWO-LEVEL MODEL OF STATIC DEFORMATION OF CONSTRUCTIONS

MADE OF REINFORCED MATERIALS

© 2015 г. О.К. Казначеева, Ю.В. Полинко

Казначеева Ольга Константиновна - канд. техн. наук, доцент, кафедра «Общеинженерные дисциплины», ЮжноРоссийский государственный политехнический университет (НПИ) имени М.И. Платова, г. Новочеркасск, Россия. E-mail: kazn_olga@mail.ru

Полинко Юлия Викторовна - аспирант, кафедра «Общеинженерные дисциплины», Южно-Российский государственный политехнический университет (НПИ) имени М.И. Платова, г. Новочеркасск, Россия. E-mail: poli-yuliya@ yandex.ru

Kaznacheeva Olga Konstantinovna - Candidate of Technical Sciences, assistant professor, department «General Engineering Disciplines», Platov South-Russian State Polytechnic University (NPI), Novocherkassk, Russia. E-mail: kazn_olga@mail.ru

Polinko Julia Viktorovna - graduate student, department of «Gneral engineering disciplines», Platov South-Russian State Polytechnic University (NPI), Novocherkassk, Russia. E-mail: poli-yuliya@yandex.ru

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

Ключевые слова: идентификация; математическая модель; композитные материалы; армированные материалы; неразрушающий контроль.

A two-level model consisting of the final element model of static deformation (the basic model) and the reduced model with explicit expressions for response functions is offered. The model represents an inhomogene-ous design, with its parameters of elasticity and rigidity representing Piecewise constant functions of coordinates, with external actions being given functions of coordinates and the response being the fields of transitions, strains and stresses. The mathematical model looks like a variation boundary-value problem with variable parametres. For practical calculations it is substituted by the discrete analogue i. e. by the final-element model. The increase of calculation profitability is obtained by introducing the reduced model in the form of explicit approximation of response function dependence on structural parametres. The reduction algorithm of the discrete model is based on the factoral computing experiment. The obtaied reduced model is used for minimization of the quality estimation criterion and for deriving interval estimations ofparametres.

Keywords: identification; mathematical model; composite materials; reinforced materials; non-destructive testing.

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

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

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

ванных теорий механики конструкций и установленных экспериментальных фактов, что исключает необходимость дополнительной настройки модели [1].

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

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

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

Натурный объект представляет собой конструкцию заданной формы и размеров, занимающий в пространстве область ОеR3. Эта область содержит заданное число непересекающихся включений - подобластей О^ с О , 01 и02 и... = 0, из материалов с различающимися физико-механическими свойствами. Взаимное расположение подобластей задано, а их размеры I(г), влияющие на интегральную жесткость, могут варьироваться, причем, не все размеры могут быть измерены без нарушения целостности изделия.

Состояние объекта характеризуется полями перемещений и(х1,х2,х3), деформаций е(х1 , Х2 , Х3 ) и напряжений с(х1, х2, х3). При этом указанные параметры состояния взаимосвязаны.

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

и еи, (1)

где и - множество всех кинематически допустимых полей перемещений.

Между перемещениями и деформациями существует взаимосвязь, точно выражаемая кинематически-

ми соотношениями: е = Уеи , где Уе - дифференциальный оператор, устанавливаемый кинематикой деформирования объекта.

Напряжения и деформации в каждой подобласти 0г- связаны определяющим соотношением материала этой подобласти. Для линейно-упругого материала определяющим соотношением является обобщенный закон Гука: с(г) = Е(г(г))е , в котором матрица упругости Е зависит от параметров d(г) - модулей упругости на растяжение-сжатие в направлениях осей анизотропии, коэффициентов Пуассона и модулей сдвига. Для физически нелинейных материалов определяющие соотношения могут быть представлены в виде

а(° = Е(0^ (г),е)е,

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

К объекту прикладываются воздействия - поле внешних сил г (х1, х2, х3). В результате возникают поля перемещений, деформаций и напряжений q={и(х1, х2, х3), е(х1, х2, х3), с(х1, х2, х3)}.

Математическая модель в форме краевой задачи, по формализации В.В. Болотина [5], отображает воздействия (поле нагрузок г) на состояние q. Оператор отображения зависит от всех структурных параметров р={1(г), d(г)}, определяющих внутреннюю геометрию армированной конструкции и характеристики упругости материалов. В лагранжевой вариационной постановке эта модель имеет вид [6]:

Н(р): г^q |ттП(р;q,u,е(u),c[е(u)]), (2)

иеи

где П - функционал потенциальной энергии, минимум которого находится при главных граничных условиях (1).

В натурном эксперименте к объекту прикладываются воздействия г и измеряются переменные состояния q. Требуется определить такие значения структурных параметров р, при которых вычисленный отклик наименее отклоняется от фактически измеренного q*:

р*:тш|р) - q * .

р

Здесь знаком || || обозначена подходящая норма (критерий качества идентификации).

На практике измерительная система позволяет измерить перемещения и деформации в местах размещения датчиков; результаты измерений обозначим через Z*. Модель (2) тоже позволяет, определив перемещения и деформации во всех точках объекта, найти их значения в местах расположения датчиков измери-

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

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

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

Модель измерительной системы должна учитывать детерминированную составляющую и случайную погрешность измерений.

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

7 = 2 + А2 = Сд + А1 ,

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

Таким образом известны:

- форма, геометрические размеры и структура объекта ^ и02 и... = 0;

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

Н(р): г ^д |ттП(р;д,п,е(п),ст[е(п)]);

п^и

- математическая модель измерительной системы

7 = Сд + А7 ;

- приложенные к натурному объекту нагрузки Г , , Х3 ) ;

- результаты натурных измерений

- критерий качества идентификации

ф(д, д*)=| |д( р) - д *.

Ограничения: двусторонние ограничения на определяемые параметры

I — i — I ршт — р — ртах .

Требуется определить:

- значения структурных параметров

p*:min| |д( р) - д *,

р

отвечающие минимуму критерия качества;

- доверительные интервалы, покрывающие с заданной вероятностью действительные значения структурных параметров с учетом случайной составляющей измерений А7 .

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

Статическое деформирование конструкции описывается на основе теорий, содержащих следующие группы уравнений [7, 8]: кинематические, в соответствии с которыми выражается взаимосвязь деформаций с полями перемещений в конструкции; уравнений равновесия, выражающих взаимосвязь между

Объект Воздействия г(хь х2, х3) -^ Результаты измерений Z *

Структурные параметры р Переменные Измерительная система *

состояния q(xi, хъ Х3) -►

К постановке задачи идентификации

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

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

В основной модели статическое деформирование описывается системой алгебраических уравнений высокого порядка:

Н(р^ = г , (3)

где Н(р) - матрица жесткости; q - вектор обобщенных перемещений; г - вектор внешних воздействий; р -определяемые параметры.

Матрица жесткости имеет следующую структуру:

H = |Фт(х) Dт ED Ф(х) dV,

(4)

H

.={Фт(х) D т Ecr D Ф(х) dV.

(5)

H 0

q = r, q,r eRn

(6)

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

Положим, что расчетная модель измерительной системы позволяет выразить детерминированную составляющую вектора измерений Z через обобщенные перемещения q: Z(p) = Cq(p).

Эта связь в дальнейшем предполагается линейной.

Критерий качества идентификации должен оценивать отклонение показаний датчиков, рассчитываемых по модели в соответствии с (6), от фактически измеренных показаний датчиков Простейший из таких критериев может быть построен как взвешенная сумма отклонений рассчитанных и измеренных величин [9]:

где интегрирование проводится по объему конструкции. Здесь Е - матрица упругости материала; Ф(х) -вектор-столбец базисных интерполяционных функций; D - дифференциальный оператор, входящий в кинематические соотношения (связь перемещений с деформациями):

5е(х) = Б 5и(х), 5и(х) = Ф(x)5q ;

5и (х) - вариация перемещений, зависящая от координат точки.

Для нелинейно-упругого материала в (4) вместо матрицы упругости Е фигурирует матрица секущих модулей Есг, зависящая от деформаций:

ф = (Z* -Cqf V(Z* -Cq) ^ min,

(7)

Матрицу коэффициентов системы (3) можно представить в виде линейной комбинации нескольких матриц, причем, коэффициентами этой комбинации будут варьируемые элементы матрицы упругости:

( V \

Здесь р, - переменные параметры, входящие в матрицу Е; Н0 - «постоянная» составляющая матрицы жесткости, рассчитываемая по формулам (4) или (5) при некоторых номинальных значениях параметров упругости; Нгр, - составляющие матрицы жесткости, пропорциональные переменным параметрам р,-.

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

где V - матрица весовых коэффициентов.

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

Минимум критерия (7) отыскивается по варьируемым параметрам р, причем, функция отклика q находится решением системы уравнений высокого порядка (3) - с постоянной матрицей коэффициентов (4), или, в случае нелинейно-упругого материала, нелинейной системы с матрицей коэффициентов (5).

Условие минимума критерия (7) приводится к

виду

(н0 - у-1но1ст уф=г - V-1н о1ст V г*.

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

Пусть для функции отклика у(х), где х - варьируемые факторы, а у - переменная состояния, построена аппроксимация

у(х) = Ль(х) + ЛХ1 (х) + Ль(х) +... + АX,(х), (8)

где х (х) - базисные функции.

После подстановки (8) в критерий (7) он принимает вид функции 5 переменных А1, ..., А 5. Если число базисных функций невелико, это делает возможным применение точного метода отыскания минимума функции. Задача заключается в построении аппроксимирующего выражения (8) и вычислении коэффициентов А,-.

Зависимость обобщенных перемещений от варьируемых параметров может быть точно представлена дробно-рациональной функцией, что позволяет приближенно представить её в виде произведения отрезков рядов Лорана [1, 10]:

q(p)=Ao + Ai

+Ac

+A

Pi P2 2

- +... +

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

(р ^ -Х2) *2... где Xi - корни характеристического уравнения

Л* (H о + Х, Н, ) = 0.

Коэффициенты дробно-рациональной аппроксимации могут быть найдены с помощью факторного вычислительного эксперимента. Особенностью этого подхода является то, что вместо непосредственного варьирования параметров необходимо варьировать базисные функции [1]. Для вычисления коэффициентов разложения необходимо перейти к новым безразмерным факторам: ^ =—1—.

р,

Вывод

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

Литература

3

Казначеева О.К., Каледин В.О. Идентификация параметров упругости и жесткости конструкций из армированных материалов: монография. Новочеркасск: Лик, 2012. 136 с.

Голованов А.И., Бережной Д.В. Метод конечных элементов в механике деформируемых твердых тел. Казань: Изд-во «ДАС», 2001. 301 с.

Ватульян А.О. Обратные задачи в механике деформируемого твердого тела. М.: Физматлит, 2007. 223 с.

4. Ватульян А.О., Нестеров С.А. Об одном подходе к восстановлению коэффициентов переноса // Изв. вузов. Сев.-Кавк. регион. Естеств. науки. 2009. № 3. С. 39 - 43.

5. Болотин В.В. Статистические методы в строительной механике: 2-е изд., перераб. и доп. М.: Стройиздат, 1965. 279 с.

6. Абовский Н.П., Андреев Н.П., Деруга А.П. Вариационные принципы теории упругости и теории оболочек. М.: Наука, 1978. 287 с.

7. Каледин В.О. Численно-аналитические модели в прочностных расчётах пространственных конструкций / НФИ КемГУ. Новокузнецк, 2000. 204 с.

8. Ломакин В.А. Теория упругости неоднородных тел. М.: Изд-во МГУ, 1976. 368 с.

9. Калашников С.Ю., Казначеева О.К., Бурцева О.А. Алгоритмы оптимального оценивания состояния и внешних воздействий наблюдаемых конструкций // Вестн. Вол-гогр. гос. архит.-строит. ун-та. Сер.: Стр-во и архит. 2011. Вып. 21(40). С. 5 - 12.

Kaledin V.O., Kaznacheeva O.K., Burtseva O.A. Stiffness

10

Parameters Identification of Nonlinear Spring Armored Beam. In: Physics and Mechanics of New Materials and Their Applications, Nova Science Publishers, NY, 2012.

References

1. Kaznacheeva O.K., Kaledin V.O. Identifikatsiya parametrov uprugosti i zhestkosti konstruktsii iz armirovannykh materialov [Identification of Elasticity and Rigidity Parameters of Structures Made of Reinforced Materials]. Novocherkassk, Lik Publ., 2012, 136 p.

2. Golovanov A.I., Berezhnoi D.V. Metod konechnykh elementov v mekhanike deformiruemykh tverdykh tel [The Finite Element Method in the Mechanics of Strained Solids]. Kazan, DAS Publ., 2001, 301 p.

3. Vatul'yan A.O. Obratnye zadachi v mekhanike deformiruemogo tverdogo tela [Inverse Problems in Solid Mechanics]. Moscow, Fizmatlit Publ., 2007, 223p.

4. Vatul'yan A.O., Nesterov S.A. Ob odnom podkhode k vosstanovleniyu koeffitsientov perenosa [An approach to the restoration of transport coefficients]. Izvestiya vuzov. Severo-Kavkazskii region. Estestvennye nauki, 2009, no. 3, pp. 39-43. [In Russ.]

5. Bolotin V.V. Statisticheskie metody v stroitel'noi mekhanike [Statistical Methods in Structural Mechanics]. Moscow, Stroiizdat Publ., 1965, 279 p.

6. Abovskii N.P., Andreev N.P., Deruga A.P. Variatsionnye printsipy teorii uprugosti i teorii obolochek [Variational Principles of Elasticity and Shell Theories]. Moscow, Nauka Publ., 1978, 287 p.

7. Kaledin V.O. Chislenno-analiticheskie modeli v prochnostnykh raschetakh prostranstvennykh konstruktsii [Numerical and Analytical Models in Strength Calculations of Spatial Structures]. Novokuznetsk, NFI KemGU, 2000, 204 p.

8. Lomakin V.A. Teoriyayprugosti neodnorodnykh tel [Theory of Inhomogeneous Bodies Elasticity]. Moscow, MSU, 1976, 368 p.

9. Kalashnikov S.Yu., Kaznacheeva O.K., Burtseva O.A. Algoritmy optimal'nogo otsenivaniya sostoyaniya i vneshnikh vozdeistvii nablyudaemykh konstruktsii [Algorithms for Optimal Estimation of the State and of the External Effects for the Observed Structures]. Vestn. Volgogr. gos. arkhit.-stroit. un-ta. Ser.: Str-vo i arkhit., 2011, vol. 21 (40), pp. 5-12. [In Russ.]

10. Kaledin V.O., Kaznacheeva O.K., Burtseva O.A. Stiffness Parameters Identification of Nonlinear Spring Armored Beam in Physics and Mechanics of New Materials and Their Applications, Nova Science Publishers, NY, 2012.

Поступила в редакцию

1 октября 2015 г.

1

1

1

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