Полный текст

Введение В последние десятилентия повысился интерес к исследованиям механических колебаний пластин в неподвижной вязкой жидкости (газе). Это обусловлено множеством практических приложений, охватывающих разнообразные области знаний, включая атомную микроскопию [1, 2], датчики и приводы головок на микромеханических генераторах [3], охлаждающие устройства [4], морскую и шельфовую технику [5, 6]. Второе направление связано с развиваемым в последнее время подходом к определению демпфирующих свойства материалов на основе исследования затухающих изгибных колебаний консольно-закрепленных плоских тест-образцов [7-9]. При этом тест-образцы находятся в воздухе, что приводит к необходимости исключения аэродинамической составляющей демпфирования с целью получения параметра внутреннего демпфирования тест-образца. Одной из основных задач в этом классе проблем является прогнозирование сил, действующих на колеблющуюся пластину со стороны жидкости (газа). Считается, что аэродинамическое взаимодействие может быть сведено к инерционному эффекту присоединенной массы и аэродинамическому демпфированию (см., например, [10-13]). Инерционный эффект приводит к снижению частоты, а аэродинамическое демпфирование - к росту декремента колебаний пластины по сравнению с ее колебаниями в вакууме. В работе [14] теоретико-экспериментальным методом построена структурная формула для определения аэродинамической составляющей демпфирования консольно-закрепленных тест-образцов прямоугольного сечения при изгибных колебаниях по низшей моде. В качестве материала тест-образцов использовался дюралюминий. У тест-образцов из такого материала внутреннее демпфирование практически не зависит от амплитуды колебаний вплоть до очень больших деформаций [15-18], что позволяет легко выделить аэродинамическую составляющую демпфирования из полного уровня демпфирования тест-образца. При этом тест-образцы находились на значительном расстоянии от окружающих предметов, чтобы они не могли повлиять на результаты эксперимента. Однако колебания конструкций возможны вблизи большеразмерных объектов, которые могут создавать экранирующий эффект, увеличивая таким образом влияние аэродинамических сил на демпфирование конструкции. С целью оценки данного эффекта актуально проведение параметрических экспериментальных исследований затухающих изгибных колебаний дюралюминиевых консольно-закрепленных тест-образцов, расположенных на различном расстоянии от абсолютно жесткого экрана. Это позволяет ставить и решать самостоятельную задачу построения методики определения аэродинамической составляющей демпфирования тест-образца (прямоугольной в плане удлиненной пластины) при наличии жесткого экрана. Однако следует заметить, что решение задачи определения сил аэродинамического сопротивления при колебаниях удлиненной пластины наталкивается на серьезные трудности (даже при отсутствии экрана), обусловленные сложностью моделирования трехмерных течений газа, вызванных колебанием пластины. Известные подходы [1, 10] к ее решению основываются на предположении, что длина пластины существенно превышает ее ширину и толщину. В этом случае на низких структурных модах колебаний длина вибрационной волны значительно больше отклонений пластины, в силу чего она может рассматриваться как локально плоская. При этом трехмерными явлениями, относящимися к течению газа вдоль длины пластины, в том числе сходом вихрей с ее торца, пренебрегают, определяя аэродинамические силы в каждом сечении пластины путем изучения плоского движения газа, вызванного гармоническими осцилляциями тонкой жесткой пластины. Однако даже в плоском приближении задача определения аэродинамических сил, действующих на гармонически колеблющуюся пластину, в полном объеме не решена. Поэтому реальным будет подход, основанный на модифицировании отмеченной в работе [14] структурной формулы с использованием экспериментальных данных о демпфирующих свойствах тест-образцов, расположенных на различных расстояниях от жесткого экрана. Для реализации данной возможности проведена модернизация разработанной ранее экспериментальной установки [7] с введением в нее абсолютно жесткого экрана, ширина которого значительно превышает ширину тест-образца (рис. 1). Рис. 1. Принципиальная схема экспериментальной установки Fig. 1. Schematic diagram of the experimental setup Установка состоит из основания 1 и силовой стойки 2, жестко соединенных между собой. На стойке неподвижно закреплена консоль 3 с захватом 4 на конце. Защемление тест-образца 5 осуществляется с помощью разнесенных жестких планок, соединенных с консолью болтовыми соединениями и исключающих поворот тест-образца в сечении заделки. На стойке установлена подвижная платформа 6 для установки лазерного датчика перемещений 7, положение которой вдоль стойки может изменяться для измерения прогиба w свободного конца тест-образца при изменении его стрелы вылета. Вблизи тест-образца расположен жесткий экран 8. Расстояние от недеформированного тест-образца до экрана может меняться. Колебания совершаются после начального отклонения свободного конца тест-образца до контакта его с жестким экраном. В установке используется триангуляционный лазерный датчик фирмы RIFTEK (RF603-X/100), обеспечивающий точность измерения амплитуды колебаний 0,01мм в цифровом формате. Разработанное математическое обеспечение позволяет осуществлять до 2000 замеров прогиба в секунду, что обеспечивает высокую точность описания экспериментальных виброграмм затухающих колебаний исследуемых тест-образцов. 1. Обработка результатов эксперимента Демпфирующие свойства тест-образца при колебаниях в воздухе определяются логарифмическим декрементом колебаний (ЛДК) d(A), зависящим от амплитуды колебаний A его свободного конца. Учитывая малость сил внутреннего трения и аэродинамических сил по сравнению с упругими силами, можно считать [8], что полный ЛДК тест-образца складывается из двух независимых частей, отвечающих механическому (внутреннее+конструкционное) dm(A) и аэродинамическому da(A) демпфированию: . (1.1) Внутреннее демпфирование обусловлено необратимым рассеянием энергии в материале [19-28], а конструкционное - потерями энергии в узле крепления тест-образца. Проведенные авторами экспериментальные исследования при различной степени зажатия тест-образцов показали, что их демпфирующие свойства остаются при этом практически неизменными. Поэтому далее считается, что параметр внутреннего демпфирования dm(A) тест-образца определяется только рассеянием энергии (внутренним трением) в материале. Обработка результатов испытаний тест-образца имеет своей целью нахождение по экспериментальной виброграмме w(t) затухающих изгибных колебаний тест-образца зависимости d(A) и низшей циклической частоты , где t - текущее время; T - период колебаний тест-образца. Обработка состоит из двух этапов. На первом этапе из записанной экспериментальной виброграммы w(t) колебаний свободного конца тест-образца выбирается рабочий диапазон изменения амплитуд [Amax; Amin] и находятся моменты времени , соответствующие нулевым прогибам w в указанном диапазоне амплитуд. По значениям t1 и tm находится осредненный полупериод колебаний , что дает возможность определить циклическую частоту . На втором этапе находятся амплитуды колебаний A1, A2,…, An в выбранном диапазоне [Amax; Amin] в области положительных значений w(t) и соответствующие им моменты времени , в результате чего получается экспериментальная дискретная зависимость ( ). Полученная зависимость аппроксимируется суммой двух экспонент: . (1.2) Параметры a1, a2, a3, a4 зависимости (1.2) определяются из условия , (1.3) где Ai - значения A(t), найденные по зависимости (1.2) в моменты времени ti. Для нахождения данных параметров можно использовать необходимое условие существования минимума целевой функции : . Однако это приводит к системе нелинейных уравнений относительно параметров a1, a2, a3, a4, для решения которой необходимо использовать итерационные методы: метод простой итерации, метод Ньютона и др. [29, 30]. При этом необходимо иметь начальные значения параметров a1, a2, a3, a4. Однако проведенные численные эксперименты показали, что выбрать данные параметры так, чтобы они обеспечивали сходимость отмеченных методов, довольно проблематично. Поэтому более подходящими являются прямые методы поиска, основанные на целенаправленном движении к минимуму целевой функции с использованием какой-либо шаговой процедуры. Из них предпочтительными будут прямые методы поиска нулевого порядка (симлекс-метод, метод Хука-Дживса, метод Розенброка) [31, 32], в которых, в отличие от градиентных методов, достаточно располагать лишь возможностью вычисления целевой функции в любой точке ее определения. Наиболее простым и удобным является метод Хука-Дживса, который легко реализуется при любой размерности пространства поиска. Для практической реализации процедуры метода Хука-Дживса предлагается использовать матрицу поиска . (1.4) Первая строка матрицы P содержит координаты a1,0, a2,0, a3,0, a4,0 базовой точки. Остальные строки содержат координаты точек, отстоящих на расстоянии шага d от базовой точки по каждому искомому параметру. Поиск осуществляется определением целевой функции по содержимому каждой строки матрицы P и выбором строки, которая дает минимальное значение . Содержимое этой строки переписывается в первую строку матрицы P и становится новой базовой точкой. Затем формируются все остальные строки матрицы P, как показано в выражении (1.4), и описанная процедура повторяется. Если текущее минимальное значение в процессе поиска соответствует первой строке матрицы P, то элементы данной строки принимаются за искомые параметры a1, a2, a3, a4 зависимости (1.2). Обработаны экспериментальные виброграммы затухающих изгибных колебаний вертикально закрепленных дюралюминиевых тест-образцов с длинами L = 200, 300 и 500 мм при расстоянии до экрана S = 70, 21, 14 и 7 мм (по два замера на каждом тест-образце). Ширина и толщина тест-образцов: b = 20 мм, h = 0,95 мм. Выбор дюралюминия продиктован тем, что данный материал обладает стабильно низкими и практически не зависящими от уровня деформации демпфирующими свойствами [15-18], что позволяет просто и надежно выделить из полного (экспериментального) ЛДК тест-образца необходимую аэродинамическую составляющую демпфирования da(A). Поиск минимума целевой функции осуществлялся методом Хука-Дживса в три этапа с уменьшением шага поиска в 10 раз после достижения очередного минимума целевой функции. Начальный шаг поиска по каждому из параметров составлял величину d = 0,001. Координаты базовой точки определялись по зависимости (1.2) при t = t1, t = tn и условиях , , что дает формулы для вычисления значений a2,0 и a1,0: . С целью получения ЛДК тест-образцов приблизительно в одинаковом диапазоне амплитуд (при различных расстояниях S до экрана) отсчет времени в аппроксимации (1.2) осуществлялся от момента, соответствующего экспериментальной амплитуде A* = 7 мм. При поиске параметров учитывалось ограничение-равенство . В табл. 1 приведены длины L тест-образцов, соответствующие им экспериментальные частоты f, расстояния S до экрана, номера замеров виброграммы колебаний w(t) и полученные из условия (1.3) параметры зависимости (1.2). На рис. 2 приведена огибающая колебаний свободного конца тест-образца длиной L = 200 мм при S = 14 мм (по первому замеру), показывающая хорошее качество поиска параметров . Таблица 1 Длины L тест-образцов, экспериментальные частоты f, расстояния S до экрана, номера замеров виброграммы w(t) и найденные параметры a1, a2, a3, a4 аппроксимации (1.2) Table 1 The length L of the test specimens, the experimental frequencies f, the distances S to the shield, the measurement numbers of the vibrogram w(t) and the found approximations (1.2) parameters a1, a2, a3, a4 L, мм f, Гц S, мм Замер a1 a2 a3 a4 200 17,670 7 1 3,51335 1,06177 3,40665 0,24345 2 3,73477 1,13263 3,25523 0,21404 14 1 3,34917 0,99389 3,66083 0,23672 2 3,60104 0,90660 3,27896 0,18927 21 1 3,08032 0,95682 3,80968 0,25265 2 4,03986 0,83692 2,88014 0,17525 70 1 3,63186 0,79979 3,28814 0,19317 2 3,12104 0,90040 3,73896 0,20828 300 8,237 7 1 3,29020 0,52347 3,06980 0,08382 2 3,58239 0,54937 3,36761 0,09576 14 1 3,76819 0,38612 3,14181 0,06693 2 3,59732 0,40716 3,36268 0,07628 21 1 3,59148 0,36400 3,26852 0,06407 2 3,71410 0,35826 3,26590 0,06789 70 1 3,49727 0,33585 3,36273 0,05893 2 3,56620 0,33794 3,35380 0,06463 500 3,148 7 1 3,17450 0,18051 3,70540 0,03102 2 2,94610 0,23521 3,88390 0,03806 14 1 3,16135 0,13636 3,71865 0,02591 2 2,90944 0,17703 3,90056 0,02932 21 1 2,90499 0,15037 3,95501 0,02550 2 2,75360 0,15857 4,11641 0,02701 70 1 2,66436 0,14355 4,27565 0,02623 2 2,67970 0,13021 4,31031 0,02937 Рис. 2. Огибающая колебаний свободного конца тест-образца длиной L = 200 мм при S = 14 мм (по первому замеру): точки - эксперимент; линия - аппроксимация (1.2) Fig. 2. The envelope of free end oscillations of a test specimen with a length L = 200 mm at S = 14 mm (according to the first measurement): points show the experiment; the line shows the approximation (1.2) Для определения ЛДК тест-образца использовалась формула [7] . После подстановки сюда аппроксимации (1.2) получается зависимость . (1.5) Выражения (1.2) и (1.5) в параметрическом виде (через время t) дают необходимую амплитудную зависимость d(A) ЛДК тест-образца. Однако это вызывает затруднения при построении усредненных по замерам зависимостей d(A) тест-образцов, обусловленные невозможностью исключения времени t из зависимостей (1.2) и (1.5) с целью получения явной зависимости d(A) на каждом замере виброграммы затухающих колебаний тест-образца. Для получения такой зависимости предлагается использовать аппроксимацию . (1.6) Коэффициенты полинома (1.6) определяются из условия . (1.7) Значения di и Ai находятся соответственно по зависимостям (1.5) и (1.2) при t = ti. Реализация условия (1.7) приводит к системе линейных алгебраических уравнений , где . Численные эксперименты показали, что наилучшее приближение к исходной зависимости (1.5) получается при использовании полинома восьмой степени (k = 8). Таким образом, можно построить аппроксимирующий степенной полином (1.6) при каждом замере виброграммы колебаний тест-образца. Усреднение полученных зависимостей достигается усреднением соответствующих коэффициентов полиномов. С целью получения зависимостей d(A) при различных замерах в одинаковых диапазонах амплитуд диапазон времени при каждом замере необходимо назначать так, чтобы выполнялись условия . (1.8) Решение уравнений (1.8) необходимо итерировать. Для этого проще всего использовать метод половинного деления с последовательным уменьшением вдвое интервала поиска до выполнения условий , (1.9) где ε - заданная точность. На рис. 3 приведены зависимости d(A) для тест-образца длиной L = 300 мм при расстоянии до экрана S = 21 мм (замер 1) в диапазоне амплитуд [7; 0,86] мм, полученные при n = 20 и ε = 1·10-6: точки - зависимость (1.5); линия - полином (1.6) при k = 8. На рис. 4 приведены исходные при замерах и усредненная по ним зависимости d(A) для тест-образца длиной L = 200 мм при S = 14 мм, полученные с использованием аппроксимации (1.6): 1 - замер 1; 2 - замер 2; 3 - усредненная зависимость. На рис. 5 приведены усредненные по замерам зависимости d(A) для тест-образца длиной L = 300 мм при расстоянии до экрана S = 7, 14, 21 и 70 мм. Видно, что по мере возрастания расстояния S ЛДК тест-образца существенно падает. Рис. 3. Зависимости d(A) для тест-образца длиной L = 300 мм при расстоянии до экрана S = 21 мм, замер 1: точки - зависимость (1.5); линия - полином (1.6) Fig. 3. Dependences of d(A) for the test specimen of length L = 300 mm at a distance to the shield S = 21 mm, measurement 1: points show the dependence (1.5); the line shows the polynomial (1.6) 3 2 1 Рис. 4. Исходные при замерах и усредненная по ним зависимости d(A) для тест-образца длиной L = 200 мм при S = 14 мм: 1 - замер 1; 2 - замер 2; 3 - усредненная зависимость Fig. 4. Initial measurement values and averaged d(A) dependencies for the test specimen of length L = 200 mm at S = 14 mm: 1 is the measurement 1; 2 is the measurement 2; 3 is the averaged dependence 1 2 3 4 Рис. 5. Усредненные по замерам зависимости d(A) для тест-образца длиной L = 300 мм при различных расстояниях S до экрана: 1 - S = 7 мм; 2 - S = 14 мм; 3 - S = 21 мм; 4 - S = 70 мм Fig. 5. Dependences of d(A) averaged over the measurements of the for the test specimen with the length L = 300 mm at different distances S to the shield: 1 - S = 7 mm; 2 - S = 14 mm; 3 - S = 21 mm; 4 - S = 70 mm 2. Построение структурной формулы для определения аэродинамической составляющей демпфирования тест-образца при наличии жесткого экрана В работе [14] при исследовании затухающих изгибных колебаний дюралюминиевых пластин теоретико-экспериментальным методом построена структурная формула для определения аэродинамической составляющей демпфирования тонкой консольно-закрепленной пластины постоянной толщины при колебаниях ее по низшей моде вдали от экранирующих объектов: (2.1) Здесь ρa, ν - соответственно плотность и кинематическая вязкость воздуха; rm - плотность материала; b - ширина пластины. Первое слагаемое в скобках отражает вклад в da вязкой составляющей аэродинамического сопротивления пластины, которая зависит только безразмерной частоты (параметра Стокса) b. Его влияние на величину da сказывается главным образом при малых амплитудах колебаний пластины. Второе слагаемое определяет вклад в da вихревой компоненты аэродинамического сопротивления, которая зависит в основном от безразмерной амплитуды κ. Для определения аэродинамической составляющей ЛДК тест-образца при наличии абсолютно жесткого экрана предлагается использовать теоретико-экспериментальный подход, состоящий в подходящем модифицировании формулы (2.1) при минимальном расхождении между расчетными и экспериментальными значениями da при нескольких амплитудах колебаний Ai ( в заданном диапазоне амплитуд [Amax; Amin]. Для учета влияния на da жесткого экрана в формулу (2.1) предлагается вести дополнительно параметры x1, x2, x3: , (2.2) определяемые из условия , (2.3) где - экспериментальная аэродинамическая составляющая демпфирования тест-образца при амплитуде колебаний Ai, получаемая вычетом из полного ЛДК параметра внутреннего демпфирования, обусловленного внутренним трением в материале тест-образца. В работе [33] показано, что в случае независимости демпфирующих свойств материала от уровня деформации параметр внутреннего демпфирования тест-образца является тождественным ЛДК dm материала. Для дюралюминия эта величина составляет значение dm = 0,0051 [17]. В табл. 2 приведены длины L тест-образцов, расстояния S до жесткого экрана и найденные из условия (2.3) параметры x1, x2, x3. Поиск минимума функции осуществлялся методом Хука-Дживса с шагом d = 0,0001 по каждому искомому параметру. Имея параметры x1, x2, x3 при четырех значениях S, можно построить кубические полиномы (при имеющихся в эксперименте длинах L тест-образцов): . (2.4) В табл. 3 приведены длины L тест-образцов, соответствующие им низшие частоты f и полученные коэффициенты ai, bi, ci (i = 0, 1, 2, 3) полиномов (2.4). Таблица 2 Длины L тест-образцов, расстояния S до экрана и найденные из условия (2.3) параметры x1, x2, x3 Table 2 The length L of the test specimens, the distance S to the shield and the parameters x1, x2, x3 found from condition (2.3) S, мм L = 200 мм L = 300 мм L = 500 мм x1 x2 x3 x1 x2 x3 x1 x2 x3 7 1,8061 1,2311 0,7104 0,8222 1,4708 0,6980 0,7062 1,2660 0,4924 14 1,7127 0,9879 0,7037 0,4616 1,0173 0,6033 0,4339 0,9036 0,4526 21 1,6791 0,8587 0,4594 0,3507 0,9257 0,6205 0,3795 0,9922 0,6931 70 1,6170 0,8277 0,6586 0,2823 0,8866 0,7321 0,4852 0,6938 0,5559 Таблица 3 Длины L тест-образцов, соответствующие им низшие частоты f и полученные коэффициенты ai, bi, ci (i = 0, 1, 2, 3) полиномов (2.4) Table 3 The length L of the test specimens, the corresponding lower frequencies f and the resulting coefficients ai, bi, ci (i=0, 1, 2, 3) of the polynomials (2.4) i L = 200 мм, f = 17,670 Гц L = 300 мм, f = 8,237 Гц L = 500 мм, f = 3,148 Гц ai bi ci ai bi ci ai bi ci 0 1,9772 1,6159 0,3776 1,5073 2,3997 0,9420 1,2632 2,2407 0,9276 1 -0,6168 -1,3280 1,5331 -2,4923 -3,4411 -0,9461 -2,0620 -3,8130 -1,9179 2 0,3900 0,6906 -1,8019 1,6298 2,4034 0,7621 1,4350 3,1573 2,0837 3 -0,0695 -0,1073 0,3962 -0,2908 -0,4411 -0,1454 -0,2598 -0,6269 -0,4475 По значениям ai, bi, ci (i = 0, 1, 2, 3) при трех частотах f построены квадратичные полиномы, определяющие зависимость ai, bi, ci от безразмерной частоты : (2.5) где Ui, Vi, Wi - векторы, содержащие коэффициенты данных полиномов (табл. 4). Таблица 4 Элементы векторов Ui, Vi, Wi (i = 0, 1, 2, 3) Table 4 Elements of vectors Ui, Vi, Wi (i = 0, 1, 2, 3) U0 1,1155´100 1,7444´10-3 1,7898´10-7 U1 -1,2898´100 -1,1502´10-2 2,7441´10-5 U2 1,0115´100 6,4248´10-3 -1,6434´10-5 U3 -1,8786´10-1 -1,0972´10-3 2,8617´10-6 V0 1,9382´100 4,5330´10-3 -1,1072´10-5 V1 -3,7736´100 -1,6968´10-3 1,4616´10-5 V2 3,5640´100 -4,5725´10-3 -3,2374´10-6 V3 -7,4384´10-1 1,4022´10-3 -1,0882´10-6 W0 8,0860´10-1 1,9483´10-3 -6,0679´10-6 W1 -2,3907´100 5,0484´10-3 6,9587´10-6 W2 2,8796´100 -9,3825´10-3 -1,1731´10-6 W3 -6,3785´10-1 2,2834´10-3 -1,8862´10-7 Подставляя полиномы (2.5) в равенства (2.4), получаем аппроксимации, определяющие зависимость параметров x1, x2, x3 структурной формулы (2.2) от относительного расстояния до экрана и безразмерной частоты b: . (2.6) На рис. 6 представлены зависимости da(A) для тест-образца длиной L = 200 мм, расположенного на расстоянии S = 14 мм от жесткого экрана: пунктирная линия - эксперимент; сплошная линия - структурная формула (2.2). Наблюдается достаточная близость расчетной зависимости к экспериментальной кривой. Аналогичная ситуация наблюдается и при всех других участвующих в экспериментах значениях L и S, что свидетельствует о возможности использования формулы (2.2) для вычисления аэродинамической составляющей демпфирования консольных тест-образцов, расположенных вблизи абсолютно жесткого экрана. Рис. 6. Зависимости da(A) для тест-образца длиной L = 200 мм, расположенного на расстоянии S = 14 мм от жесткого экрана: пунктирная линия - эксперимент; сплошная линия - структурная формула (2.2) Fig. 6. da(A) dependences for the test specimen of length L = 200 mm located at a distance S = 14 mm from the rigid shield: the dashed line shows the experiment; the solid line shows the structural formula (2.2) Выводы 1. Разработана простая и удобная в практическом применении методика обработки виброграммы затухающих изгибных колебаний консольно-закрепленных тест-образцов для определения экспериментальной низшей частоты f и амплитудной зависимости ЛДК d тест-образца. Проведены численные эксперименты, подтверждающие достоверность разработанной методики. 2. Проведены экспериментальные параметрические исследования, показывающие значительное увеличение демпфирующих свойств тест-образца при уменьшении расстояния S до жесткого экрана. 3. Разработан теоретико-экспериментальный метод построения структурной формулы для вычисления аэродинамической составляющей демпфирования da тест-образца в зависимости от его низшей частоты колебаний и расстояния до жесткого экрана, основанный на модификации формулы, полученной ранее для тонкой пластины постоянной толщины при отсутствии экрана, путем введения в нее трех дополнительных параметров, определяемых из условия минимума квадратичной невязки между экспериментальными и расчетными значениями аэродинамической составляющей демпфирования тест-образца. Получены соотношения, определяющие зависимости отмеченных параметров от относительного расстояния q = S/b до экрана и безразмерной частоты колебаний тест-образца. Проведены численные эксперименты, подтверждающие достоверность разработанного метода.

Об авторах

В Н Паймушин

Казанский национальный исследовательский технический университет им. А.Н. Туполева; Казанский федеральный университет

В А Фирсов

Казанский национальный исследовательский технический университет им. А.Н. Туполева

Р К Газизуллин

Казанский национальный исследовательский технический университет им. А.Н. Туполева

В М Шишкин

Вятский государственный университет

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

Статистика

Просмотры

Аннотация - 240

PDF (Russian) - 78