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

Размещено на http://www.allbest.ru/

Кубанский государственный технологический университет

ИССЛЕДОВАНИЕ, ОЦЕНКА И ВЫБОР ПАРАМЕТРОВ ТЕХНИЧЕСКОГО СОСТОЯНИЯ ДВИГАТЕЛЕЙ ПРИ ИХ КОНТРОЛЕ И ПРОГНОЗИРОВАНИИ

Шевцов Юрий Дмитриевич д.т.н., профессор

Атрощенко Валерий Александрович д.т.н., профессор

Дудник Людмила Николаевна к.т.н.

Горохов Дмитрий Александрович к.т.н., доцент

Краснодар, Россия



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

Ключевые слова: ТЕХНИЧЕСКОЕ СОСТОЯНИЕ ДВИГАТЕЛЯ, ВЫБОР ПАРАМЕТРОВ, ПРОГНОЗИРОВАНИЕ, ФАКТОРНЫЙ АНАЛИЗ

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

Объектом исследования является двигатель, техническое состояние которого можно описать совокупностью параметров, представляющих собой вектор технического состояния , где n - число параметров системы.

Целью данной работы является повышение эксплуатационных характеристик двигателей энергетических установок (ЭУ) за счет получения достоверной информации для контроля и прогнозирования их технического состояния. Для реализации цели возникает необходимость решения следующих задач:

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

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

Исследования и оценка параметров технического состояния двигателей ДЭС

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

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

Процентная концентрация продуктов износа в масле, также как и других нерастворимых в нем загрязняющих примесей, по времени работы t может быть определена выражением (1) [1,3]:

, (1)

где

X₀ - начальная концентрация продуктов износа в масле при t=0;

e - основание натурального логарифма;

Q - скорость расхода масла через фильтр;

- коэффициент отсева фильтра;

Q_у - скорость расхода масла на угар;

V - количество масла в двигателе; - скорость поступления продуктов износа в масло;

- удельный вес продуктов износа.

1 и 3 - дизель 2Ч 10,5/13 без средств очистки с полнопоточным фильтром грубой очистки; 2 и 4 - дизель Ч8,5/11 без средств очистки и с полнопоточным фильтром тонкой очистки; 5 - дизель 6ЧН 12/14 с полнопоточным фильтром тонкой очистки

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

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

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

Методика определения необходимых зависимостей в статье рассматривается на примере гидравлической модели масляного фильтра, полученной на основе проведенных экспериментальных исследований в системе смазки шести цилиндрового дизельного двигателя 4ч14/16 с сетчатым масляным фильтром [2].

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

(2)

Модуль и аргумент выражения (2) определяют амплитудную и фазовую частотные характеристики масляного фильтра.

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

а) активная составляющая гидравлического сопротивления фильтра;

б) инерционная составляющая гидравлического сопротивления фильтра;

1-промежуточный закон; 2-закон постепенного закупоривания пор; 3-закон полного закупоривания пор; 4-закон с образованием осадка

Рисунок 2 Изменение активной и инерционной составляющих гидравлического сопротивления фильтра в соответствии с законами фильтрования

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

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

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

Рисунок 3 Изменение АЧХ фильтра во всем диапазоне времени эксплуатации до технического обслуживания

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

На рисунке 4 показано изменение амплитуд А₁(щ) ;А₂(щ); А₃(щ), а также скоростей их изменения V₁(щ); V₂(щ) и V₃(щ) на разных резонансных частотах в зависимости от времени их эксплуатации

а) б)

а) изменение амплитуд А₁(щ) А₂(щ) А₃(щ) и б) скоростей V₁(щ) V₂(щ) V₃(щ) на разных резонансных частотах от времени эксплуатации

Рисунок 4 Параметры амплитудной частотной характеристики масляного фильтра

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

На рисунке 5 показаны графики изменения диагностических параметров =, R_ф=, = и ; от времени эксплуатации.



Рисунок 5 Изменение характеристик фильтра от длительности эксплуатации

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

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

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

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

Рассмотрим применение методики выбора совокупности параметров, описывающих техническое состояние двигателей путем оценки их приоритетности, исключения малозначащих и малоинформативных параметров. Методология выбора параметров для контроля и прогнозирования основывается на использовании методов факторного анализа, а именно метода главных факторов [5].

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

R_ф - активная составляющая гидравлического сопротивления фильтра - y₁; - инерционная составляющая гидравлического сопротивления фильтра - y₂; X - концентрация примесей в масле - y₃; W - средняя производительность фильтра - y₄; А₁(щ) А₂(щ) А₃ (щ) - амплитуды амплитудно-частотных характеристик фильтра на различных резонансных частотах- y₅, y₆, y₇; V₁(щ) V₂(щ) V₃(щ) - скорость изменения амплитуд амплитудно-частотной характеристики фильтра на различных резонансных частотах - y₈, y₉, y₁₀.

Измерения параметров двигателей проводились через равные промежутки в моменты времени от t = 0 до t = 270 часов. В работе использованы следующие условные обозначения: j - номер контролируемого двигателя (j =1,…, m); i - номер параметра контроля (i =1, …., n), k - номер измерения; s - количество измерений.

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

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

Построение корреляционной матрицы осуществляется на основе расчета коэффициентов парных корреляций параметров с использованием следующих формул [6,7]:

1). Среднее значение параметров по данным их измерений при каждом замере (табл. 1):

, (3)

где i - номер параметра, j - номер объекта; k - номер измерения, m - количество контролируемых объектов.

Таблица 1

Средние значения параметров при каждом замере

Номер измере-ния	y1	y2	y3	y4	y5	y6	y7	y8	y9	y10
k=1	0,103	0,798	0,161	0,705	3,505	1,095	1,94	0,004	0,002	0,003
k=2	0,125	1,028	0,858	4,897	2,355	0,597	1,24	0,038	0,016	0,023
k=3	0,135	1,572	1,727	5,295	1,639	0,4	0,80	0,005	0,001	0,005
k=4	0,225	2,553	1,751	6,295	1,476	0,267	0,55	0,002	0,002	0,0003
k=5	0,425	4,519	1,256	9,701	1,338	0,247	0,465	0,0023	0,0003	0,0003
k=6	0,978	7,429	0,97	17,50	1,204	0,226	0,501	0,0023	0,0003	0,0002
общее среднее значен-ие	0,332	2,983	1,121	7,399	1,9195	0,472	0,919	0,0084	0,0035	0,0049

2). Произведем усреднение параметров по числу измерений s:



, (4)

где s - количество измерений.

3). Отклонения от общего среднего при каждом замере рассчитаем по формуле:

5)

4). Коэффициенты ковариации определяются по формулам:

(6)

(7)

(8)

5). На основании ковариационной матрицы строится корреляционная матрица параметров . Коэффициенты парных корреляций переменных рассчитываются по следующей формуле:

(9)

Корреляционная матрица параметров



На следующем этапе осуществляют построение редуцированной корреляционной матрицы , где i =1,…, n, k = 1,…, n с общностями на главной диагонали. Оценка общностей производится на основе вычисления максимального элемента столбца (строки) и простановки его на главной диагонали матрицы корреляций.

Редуцированная корреляционная матрица ^:

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

, (10)

где б_il - собственный вектор; л_l - собственные значения редуцированной корреляционной матрицы; i, l = 1, 2, …., n.

Для определения собственных значений л_l и собственных векторов б_il редуцированной корреляционной матрицы используем итерационный метод. Для этого необходимо вычислить суммы элементов строк матрицы , разделить каждую полученную сумму на максимальное значение суммы и получить элементы вектора б⁽¹⁾. Вектор б⁽¹⁾ перемножают с матрицей R по формуле:

(11)

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

(12)

Повторяют процедуру (11) с вектором .

(13)

где верхние индексы в скобках означают шаг итерации.

Формула (13) является общей для k шагов итерации. Процесс повторяется до тех пор, пока не добиваются сходимости к первому собственному значению R: и соответствующему первому собственному вектору . Итерационный процесс заканчивается, когда и с достаточной точностью совпадают друг с другом.

Для первого собственного вектора =( , …,), были найдены значения элементов = (-0,852; -0,797; -0,607; -0,780; 0,996; 1; 0,548; 0,572; 0,593), а также первое собственное значение л₁ = max(),
л₁ = 5,421182.

Вычислим нагрузки первого фактора a_i1 матрицы факторного отображения А=(a_il) используя формулу (10). Результаты вычислений представлены в таблице 2. Для контроля проверяем =5,421, т.е. первое собственное число матрицы равно сумме квадратов нагрузок первого фактора.

Таблица 2

Результаты вычисления первого собственного вектора , вектора и нагрузок первого фактора a_i1

=первый собственный вектор
-0,852537	-4,621730	-0,794225
-0,797670	-4,324284	-0,743111
-0,607591	-3,29387	-0,56603
-0,780137	-4,229263	-0,726777
0,996425	5,401810	0,928271
0,965319	5,233182	0,899293
1	5,421182	0,931601
0,548372	2,972737	0,510864
0,572631	3,10425	0,533464
0,593504	3,217407	0,552909

Определим матрицу воспроизведенных корреляций R⁺ с учетом только первого фактора и остаточной матрицы корреляций R₁. Матрица воспроизведенных корреляций определяется по формуле:

, (14)

где - вектор столбец нагрузок первого фактора, - вектор строка нагрузок первого фактора.

Определяем остаточную матрицу корреляций R₁, которая остается после выделения первого фактора по следующей формуле:

, (15)

где - матрица воспроизведенных корреляций ; а - редуцированная корреляционная матрица.

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

В результате вычислений получаем элементы второго собственного вектора = (-0,0472; -0,1467; 0,4293; 0,2631; -0,4622; -0,5277; -0,4165; 1; 0,9895; 0,9777) и второе собственное число матрицы : = 2,686.

Факторные нагрузки a_i2 второго фактора находим по формуле (10), и получаем следующие значения: a_i2 = -0,393; -0,122; 0,357; 0,218; -0,334; -0,3392; -0,340; 0,832; 0,823; 0,813.

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

В результате расчетов выделено три общих фактора и получена матрица факторного отображения A=(a_il), где каждый фактор характеризуется столбцом, каждая переменная - строкой матрицы А:

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

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

, (16)

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

Определим дисперсии каждого выделенного общего фактора, для чего рассчитаем квадраты факторных нагрузок по матрице факторного отображения A=(a_il). Квадраты факторных нагрузок показывают, какая часть дисперсии i-го показателя описывается k-м фактором, а сумма квадратов факторных нагрузок показывает вклад всех факторов в дисперсию каждого показателя. Характерность рассчитывают по формуле:

. (17)

Результаты вычислений приведены в таблице 3.

Таблица 3

Расчет квадратов факторных нагрузок a_il², общностей и характерности

Квадраты факторных нагрузок ail2	Общность	Характерность
F1()	F2()	F3()
0,6307936	0,001544694	0,358866962	0,991205286	0,008794714
0,5522148	0,014909009	0,131711195	0,698835052	0,301164948
0,320394	0,127663267	0,403032082	0,851089328	0,148910672
0,5282054	0,047957788	0,066063295	0,642226478	0,357773522
0,8616888	0,112020674	0,010571341	0,984280798	0,015719202
0,8087281	0,115037348	0,010860317	0,934625781	0,065374219
0,867882	0,116057601	0,011236017	0,99517521	0,00482479
0,2609828	0,692597061	0,040543603	0,994123468	0,005876532
0,2845844	0,678173297	0,024484171	0,987241855	0,012758145
0,3057093	0,662119633	0,017865187	0,98569413	0,01430587

Рассчитаем дисперсии каждого выделенного фактора как сумму квадратов факторных нагрузок:

Дисперсия первого фактора D₁= 5,421; дисперсия второго фактора D₂ = 2,56; дисперсия третьего фактора D₃ = 1,075; суммарная общность составляет 9,0645; а суммарная характерная дисперсия - 0,93.

В результате расчетов видно, что на три первых фактора приходится 90,64% от полной дисперсии, и только 9,3% приходится на суммарную характерную дисперсию. Произведем оценку значимости параметров y_i в зависимости от факторных нагрузок в первый, второй и третий общий фактор. Наибольший вклад в первый общий фактор вносит параметр Y₇ - амплитуда амплитудно-частотная характеристика фильтра А₃(щ), во второй фактор Y₈ -скорость изменения амплитудно-частотной характеристики фильтра V₁(щ) и в третий фактор параметр Y₁ - активная составляющая гидравлического сопротивления фильтра R_ф.

На основе проведенного анализа и расчетов можно заключить, что из всей совокупности характеристик достаточно выбрать параметры y_7, y_8, y₁ для прогнозирования технического состояния двигателей ДЭС. Выбранными являются параметры АЧХ, разработанной авторами модели фильтра, а именно, амплитуды А₃(щ),скорости их изменения V₁(щ) и гидравлическое сопротивления R_ф, по которым оценивают критическую степень загрязнения фильтра, определяющую необходимость проведения технического обслуживания двигателя. Остальные параметры можно не учитывать, так как они коррелированны между собой и вносят меньший вклад в общие факторы.



Заключение

технический двигатель смазка факторный

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

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



Список литературы

1. Атрощенко В.А. и др. Технические возможности повышения ресурса автономных электростанций энергетических систем: Монография [Текст] / В.А. Атрощенко, Ю.Д. Шевцов, П.В. Яцынин и др. Педько. Краснодар: Издательский Дом-Юг, 2010. 192 с.

2. Атрощенко В.А., Шевцов Ю.Д., Кобзева С.А. и др. Определение периодичности технического обслуживания двигателей внутреннего сгорания по параметрам частотной модели фильтра. Научный журнал КубГАУ [Электронный ресурс]. Краснодар: КубГАУ, 2013. №10(094). С. 704 - 717. IDA.

3. Григорьев М.А. Очистка масел двигателей внутреннего сгорания. М.: Машиностроение, 1983. 148 с.

4. Жужиков В.А. Фильтрование: Теория и практика разделения суспензий. 4-е изд., перераб. и доп. М.: Химия, 1980. 400 с.

5. Иберла К. Факторный анализ. М.: Статистика, 1980. 398 с.

6. Петриченко Г.С., Дудник Л.Н. Выбор параметров контроля технического состояния для цифровых блоков корпоративной сети на основе использования методов факторного анализа // Автоматизация и современные технологии. 2010. № 2. С. 16 - 21.

7. Петриченко Г.С., Нарыжная Н.Ю., Дудник Л.Н. Диагностика и прогнозирование технического состояния компьютерной сети. Краснодар: Издательский Дом - ЮГ, 2010. 188 с.

Размещено на Allbest.ru

...