Численное моделирование вынужденных колебаний стержня в пространственном потоке

Вестник ПНИПУ. Аэрокосмическая техника. 2014. № 39

УДК 534-143:532.527

А.В. Ермаков, И.К. Марчевский, Г.А. Щеглов

Московский государственный технический университет им. Н.Э. Баумана, Москва, Россия

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ВЫНУЖДЕННЫХ КОЛЕБАНИЙ СТЕРЖНЯ В ПРОСТРАНСТВЕННОМ ПОТОКЕ

Исследуется актуальная для многих технических приложений задача прямого численного моделирования вынужденных колебаний стержня в потоке, вызываемых его нестационарным пространственным вихревым обтеканием. Описан оригинальный алгоритм численного моделирования, в котором динамика стержня рассчитывается методом разложения по собственным формам, а процесс вихреобразования моделируется бессеточным лагранжевым методом вихревых элементов. Собственные частоты и формы колебаний определяются методом конечных элементов в коммерческом пакете MSC.Nastran. В методе вихревых элементов использована гипотеза потока завихренности, а в качестве вихревого элемента выбрана новая модель симметричного вортона-отрезка. Гидродинамические нагрузки рассчитываются с использованием аналога интеграла Коши - Лагранжа. Программная реализация метода вихревых элементов использует технологию распараллеливания MPI. Представлены результаты тестирования алгоритма на модельной задаче. Динамика стержня описывается шестью первыми тонами колебаний; вихревой след моделируется при помощи нескольких десятков тысяч вортонов. В результате расчета за стержнем формируется пространственная вихревая дорожка Кармана, частота схода вихрей близка к наблюдаемой в экспериментах. Исследовано пространственное движение стержня. Получены траектории движения сечений стержня и спектр их перемещений, в котором преобладает частота схода вихрей Кармана и вторая собственная частота колебаний. Установлен эффект снижения лобового сопротивления упругого стержня на 15-20 % по сравнению с абсолютно жестким. Разработанный алгоритм позволяет исследовать взаимодействие элементов конструкций с набегающим потоком, при этом форма и упругая модель могут быть достаточно сложными. По сравнению с известными сеточными методами предложенный подход позволяет существенно сократить время проведения расчетов.

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

A.V. Ermakov, I.K. Marchevskiy, G.A. Shcheglov

Bauman Moscow State Technical University, Moscow, Russian Federation

NUMERICAL SIMULATION OF FORCED OSCILLATIONS OF THE ROD IN THE SPATIAL FLOW

A relevant for number of technological applications problem of the direct numerical simulation of forced rod oscillations in the flow is investigated. The oscillations are caused by unsteady three-dimensional vortex flow around the rod. An original algorithm for the numerical simulation, in which the rod dynamics is simulated by using normal modes analysis method, while the vortex generation process

is simulated by using some modification of meshless Lagrangian vortex elements method. Eigenfre-quencies and eigenmodes are computed by using the finite element method in the commercial software package MSC.Nastran. In the vortex element method the hypothesis of vorticity flux is taken into account, and the original mathematical model of symmetrical vortex fragmenton is chosen as vortex element model. Hydrodynamic loads are computed using Cauchy - Lagrange integral analog. The developed vortex element method-based software uses MPI-parallelism. The results of the algorithm verification are presented for the particular model problem. The rod dynamics is described by six lowest eigenmodes; vortex wake consists of few tens of thousands of vortex fragmentons. Von-Karman-type spatial vortex street is simulated behind the rod and the vortex shedding frequency is close to that observed in the experiments. The spatial motion of the rod is investigated. The trajectories of the rod cross-sections motion are obtained as well as their movements' frequency spectra. The frequency of von-Karman vortex shedding and the second eigenfrequency dominate in the spectra. The phenomenon of the 15-20 % drag reduction for the elastic rod in comparison with the absolutely rigid one is pointed out. The developed algorithm enables to investigate the interaction of structural elements with the flow, and the shape and the elastic model can be quite complex. In comparison with well-known mesh methods the developed approach can significantly reduce the time of calculation.

Keywords: hydroelasticity, oscillations eigenmodes, finite element method, vortex element method, vortisity flux, symmetrical vortex fragmenton, parallel computing, von-Karman vortex street, oscillation spectrum, drag force.

Введение

Стержень, взаимодействующий с потоком несжимаемой среды, является распространенной расчетной схемой, применяемой при математическом моделировании широкого класса технических систем. В случае протяженного жесткого стержня постоянного кругового сечения, помещенного в поперечный поток, гидродинамические нагрузки определяются сходом вихрей Кармана и близки к моногармоническим с безразмерной частотой Sh ~ 0,2 [1]. Однако если рассматриваемый стержень является упругим, для правильного описания его вынужденных колебаний в потоке требуется решение сопряженной задачи гидроупругости. Данная задача рассматривалась многими авторами, однако на сегодняшний день прямое численное моделирование пространственного обтекания податливого стержня с использованием сеточных методов связано со значительными затратами вычислительных ресурсов.

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

Отличием предлагаемой модификации вихревого метода от известных панельных методов и метода дискретных вихрей [2] является

отсутствие необходимости явного задания линии отрыва потока от поверхности обтекаемого тела за счет использования гипотезы потока завихренности [3].

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

Рассматривается прямолинейный стержень длиной Ь, имеющий постоянное поперечное сечение в форме кольца радиусом Я с толщиной стенки И. Стержень имеет погонную массу у, материал стержня изотропный с модулем упругости Е и коэффициентом Пуассона V. Концы стержня неподвижно закреплены и заглушены. Стержень помещен в поток несжимаемой среды постоянной плотности р, имеющей

скорость Ух (рис. 1).

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

Математическая модель сопряженной задачи гидроупругости имеет вид

Математическая модель

Рис. 1. Расчетная схема стержня в поперечном потоке

У-У = 0,

^+(V .у)у=^у2 V-у (р \

д р I р)

где и - вектор перемещения профиля в сечении стержня; V - поле скоростей жидкости; р - давление. Граничные условия соответствуют заделке концов стержня, прилипанию жидкости к обтекаемой поверхности гк = г0 к + и( х,^) и затуханию возмущений на бесконечности:

I =ди

и х=-ШЛп = Эх

= 0, У(Гк= Гк , V ^ V., р ^ рх .

кк

х=-Ь/2, Ы2

Начальные условия соответствуют равновесному положению стержня в покоящейся среде.

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

Численный метод

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

Поверхность стержня заменяется системой многоугольных (чаще всего - четырехугольных) панелей. Деформация стержня описывается вектором перемещений вершин панелей . Для каждой панели заданы контрольная точка к 0, расположенная в центре панели, и вектор единичной внешней нормали п0. Вихревая рамка располагается на малом расстоянии в над панелью и состоит из т1 вортонов-отрезков, характеризуемых положениями маркеров г/ и вихревых векторов И ^,

^ = 1, ..., то.

г

Рис. 2. Расчетная схема вортонной рамки над панелью обтекаемой поверхности

Циркуляции вихревых рамок находятся из условия равенства нормальных компонент вектора скорости среды и обтекаемой поверхности в контрольных точках, аналогично методу дискретных вихрей с замкнутыми вихревыми рамками [5]. Далее в соответствии с моделью потока завихренности [6] весь вихревой слой на поверхности тела становится свободным и формирует вихревой след. В расчетах этот процесс моделируется разделением вортонных рамок на отдельные ворто-ны-отрезки, маркеры которых движутся по траекториям жидких частиц, а вихревые векторы изменяются в соответствии с уравнениями гидродинамики, как описано в работе [4]. Интенсивности вортонов-отрезков определяются в момент разделения вортонных рамок и далее остаются постоянными.

Для расчета поля скоростей среды по известным положениям и интенсивностям вортонов используется закон Био - Савара; для восстановления давления используется аналог интеграла Бернулли и Ко-ши - Лагранжа [7]. Численное интегрирование уравнений движения вортонов-отрезков производится явным методом Эйлера первого порядка точности с постоянным шагом At.

Алгоритм метода вихревых элементов допускает эффективное распараллеливание: в разработанном программном комплексе, основанном на использовании технологии MPI, доля параллельного кода доведена до 98 %, что позволяет производить расчеты на вычислительных кластерах, содержащих десятки вычислительных ядер [8-10].

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

лов сетки на поверхности } связан с вектором главных координат (ф} при помощи матрицы собственных форм: = [ А](ф}.

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

{ф}+ад^ (ф }+[ю2 ^ (ф)=(/д к

где 8 - коэффициент модального демпфирования; [ю]^ - диагональная матрица собственных частот; (/н} - вектор нагрузки, приведенный к главным координатам.

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

ф* (О=^г+

Л , ехР(-пюкДО

ю

ю

2л/т2

х{ [ю2л/Г

п ф* ^ г) - Л

С08(ю*л/1 - п2 Дt) +

+

Юкфк ^ г ) + п ( ф* ^ г ) - !к )] 1 - п2 Дt) }

ф к

ехр

п

{

п 2 фк ^ г )

С0Б(юк\/ 1 - п2 Дt) +

+

Л

юк

-

л/1 - п2 Дt) I.

Собственные частоты и формы колебаний определяются методом конечных элементов (рис. 3). В препроцессоре Ра1хап на поверхности стержня строится сетка из четырехугольных элементов, используемая для конечно-элементного расчета. Координаты узлов конечных элементов используются для задания описанных выше панелей на профиле в гидродинамическом расчете. Матрицы [ю]^ и [А] находятся

с помощью решателя БОЫОЗ пакета М8С.Кав1хап.

Тестирование алгоритмов производилось путем сравнения с известными экспериментальными результатами М.Б. Гафурова, М.А. Иль-гамова [11]. Подробное описание проведенных расчетов приведено в работе [12].

Рис. 3. Блок-схема процесса подготовки исходных данных

Вычислительный эксперимент

В качестве модельной задачи рассматривался стержень с параметрами Ь = 1,1 м, Я = 0,025 м, И = 0,001 м, изготовленный из изотропного материала (Е = 2-109 Па, у = 0,3, у = 1,1 кг/м ). Коэффициент модального демпфирования принимался 8 = 0,05. Конечно-элементная модель содержала 1792 элемента типа 0ЦАО4, что в методе вихревых элементов соответствовало 16 панелям в окружном направлении, 110 панелям вдоль образующей и 32 панелям на торцах. Скорость набегающего потока была принята равной 1 м/с.

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

ю12 = 26,6 Гц, ю34 = 70,8 Гц, ю56 = 132,7 Гц.

Рис. 4. Используемые собственные формы колебаний

В расчете проводилось моделирование переходного режима от момента времени ^ = 0 до ^ = 1с. Шаг интегрирования уравнений был выбран At = 0,001с. Вихревой след за стержнем в момент окончания расчета представлен на рис. 5, где точками обозначены положения маркеров вихревых элементов.

Рис. 5. Вихревой след за стержнем в момент времени t = 1 с

На рис. 5 заметны характерные «вихревые сгустки», соответствующие крупным вихрям Кармана. В среднем сечении след за стержнем имеет характерный вид вихревой дорожки (рис. 6).

А

Рис. 6. Дорожка Кармана за стержнем в сечении 0x2

В ходе численного моделирования установлено, что под действием набегающего потока стержень в плоскости Оху изгибается в основном по первой собственной форме (рис. 7), тогда как в плоскости 0x2 проявляется вторая собственная форма. В результате движение сечений имеет сложную траекторию (рис. 8).

Рис. 7. Деформированная форма стержня в момент времени t = 0,1 с

Рис. 8. Траектории центров сечений стержня при 0,8 < t < 1,0, соответствующих пучностям 1-й и 2-й собственных форм

На рис. 9 в качестве примера приведены зависимости от времени перемещений среднего сечения стержня (х = 0) в направлении осей Оу и 02.

Анализ поперечных перемещений сечений показал, что в направлении потока после окончания переходного режима (при t > 0,3 с) возникает статический прогиб, форма которого близка к первой форме колебаний, и вынужденные колебания. В спектре вынужденных колебаний преобладают частоты, равные 66 и 3,3 Гц, близкие ко второй собственной частоте и частоте схода вихрей Кармана соответственно.

Рис. 9. Зависимости перемещений среднего сечения стержня

На рис. 10 показан фрагмент спектра перемещений среднего сечения в окрестности частоты схода вихрей Кармана, в котором частоты

2 Я

приведены к числам Струхаля: БЬ :

0,2 (частота схода вихрей) /

0 0,2 0,4 0,6 0,8

Рис. 10. Фрагмент спектра перемещений центра среднего сечения стержня в направлении оси 02

Также установлено, что вынужденные колебания стержня приводят к некоторому снижению лобового сопротивления. На рис. 11 показаны зависимости силы лобового сопротивления для жесткой (rigid) и упругой (elastic) моделей.

0 200 400 600 800 1000 0 200 400 600

Рис. 11. Лобовое сопротивление упругого и абсолютно жесткого стержней

Среднее значение силы лобового сопротивления, показанное на рис. 11 сплошной линией, снижается при этом примерно на 15-20 %.

Отметим, что среднее время проведения расчета на четырехъя-дерных персональных ЭВМ не превышало суток. Объем необходимой оперативной памяти не превышал 10 Мб; суммарный объем дискового пространства при сохранении результатов на каждом шаге расчета не превышал 100 Мб. При проведении аналогичных расчетов сеточными методами затраты вычислительных ресурсов, как правило, получаются на порядок большими.

Заключение

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

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

Авторы благодарят Центр компетенций компании MSC Software приМГТУ им. Н.Э. Баумана. Работа выполнена при частичной финансовой поддержке гранта Президента РФ (проект МК-3705.2014.8).

Библиографический список

1. Жукаускас А., Улинскас Р., Катинас В. Гидродинамика и вибрации обтекаемых пучков труб. - Вильнюс: Мокслас, 1984. - 312 с.

2. Нелинейная теория крыла и ее приложения / Т.О. Аубакиров, С.М. Белоцерковский, А.И. Желанников, М.И. Ништ. - Алматы: Гы-лым, 1997. - 448 с.

3. Cottet G.-H., Koumoutsakos P. Vortex Methods: Theory and Practice. - Cambridge University Press, 2008. - 327 p.

4. Марчевский И.К., Щеглов Г.А. Модель симметричного ворто-на-отрезка для численного моделирования пространственных течений идеальной несжимаемой среды // Вестник МГТУ им. Н.Э. Баумана. Естественные науки. - 2008. - № 4. - С. 62-71.

5. Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент (в математической физике, аэродинамике, теории упругости и дифракции волн). - М.: Янус, 1995. - 521 с.

6. Lighthill M.J. Introduction. Boundary layer theory // Laminar Boundary Layers / Ed. L. Rosenhead. - Dover, Mineola, 1963. - P. 46-113.

7. Андронов П.Р., Гувернюк С.В., Дынникова Г.Я. Вихревые методы расчета нестационарных гидродинамических нагрузок. - М.: Изд-во МГУ, 2006. - 184 с.

8. Марчевский И.К., Щеглов Г. А. Применение параллельных алгоритмов при решении задач гидродинамики методом вихревых элементов // Вычислительные технологии и программирование. - 2010. -Т. 11. - С. 105-110.

9. Марчевский И.К., Щеглов Г.А. Параллельная реализация метода вихревых элементов с использованием модели симметричного вортона-отрезка // ПаВТ'2013: тр. междунар. науч. конф. - Челябинск, 2013. -С. 203-214.

10. Учебно-экспериментальный вычислительный кластер. Ч. 2. Примеры решения задач / В.В. Лукин, И.К. Марчевский, В.С. Морева, А.Ю. Попов, К.Л. Шаповалов, Г.А. Щеглов // Вестник МГТУ им. Н.Э. Баумана. Естественные науки. - 2012. - № 4. - С. 82-102.

11. Гафуров М.Б., Ильгамов М.А. Изгиб цилиндрической оболочки конечной длины при ее поперечном обтекании жидкостью // Прикладная механика. - 1978. - Т. 14, № 3. - С. 60-67.

12. Ермаков А. В., Щеглов Г.А. Моделирование методом вихревых элементов динамики цилиндрической оболочки в пространственном потоке жидкости // Известия вузов. Машиностроение. - 2014. -№ 3. - С. 35-41.

References

1. Zhukauskas А., Ulinskas R., Katinas V. Gidrodinamika i vibratsii obtekaemykh puchkov trub [Hydrodynamics and vibration of the streamlined bundles]. Vilnyus: Mokslas, 1984. 312 p.

2. Aubakirov T.O., Belotserkovskiy S.M., Zhelannikov А.Г, Nisht M.I. Nelineynaya teoriya kryla i ee prilozheniya [Nonlinear wing theory and its applications]. Almaty: Gylym, 1997. 448 p.

3. Cottet G.-H., Koumoutsakos P. Vortex Methods: Theory and Practice. Cambridge University Press, 2008. 327 p.

4. Marchevskiy I.K., Shheglov GA. Model simmetrichnogo vortona-otrezka dlya chislennogo modelirovaniya prostranstvennykh techeniy ideal-noy neszhimaemoy sredy [Symmetric vortex fragmenton mathematic model for the numerical simulation of three-dimensional inviscid incompressible

flows]. VestnikMoskovskogo gosudarstvennogo tekhnicheskogo universiteta imeni N.E. Baumana. Estestvennye nauki, 2008, no. 4, pp. 62-71.

5. Lifanov I.K. Metod singulyarnykh integralnykh uravneniy i chislennyy ehksperiment (v matematicheskoy fizike, aerodinamike, teorii uprugosti i difraktsii voln) [Singular integral equations method and numerical experiment (in mathematical physics, aerodynamics, the theory of elasticity and wave diffraction)]. Moscow: Yanus, 1995. 521 p.

6. Lighthill M.J. Introduction. Boundary layer theory. Laminar Boundary Layers. Ed. L. Rosenhead. Dover, Mineola, 1963. P. 46-113.

7. Andronov P.R., Guvernyuk S.V., Dynnikova G.Ya. Vikhrevye me-tody rascheta nestatsionarnykh gidrodinamicheskikh nagruzok [Vortex methods for unsteady hydrodynamic loads computation]. Moskovskiy gosu-darstvennyy universitet, 2006. 184 p.

8. Marchevskiy I.K., Shheglov G.A. Primenenie parallelnykh algorit-mov pri reshenii zadach gidrodinamiki metodom vikhrevykh elementov [Application of parallel algorithms for solving CFD-problems using vortex element method]. Vychislitel'nye tekhnologii i programmirovanie, 2010, vol. 11, pp. 105-110.

9. Marchevskiy I.K., Shheglov G.A. Parallelnaya realizatsiya metoda vikhrevykh elementov s ispolzovaniem modeli simmetrichnogo vortona-otrezka [Parallel implementation of the vortex element method using the symmetric vortex fragmenton mathematic model]. Trudy mezhdunarodnoy nauchnoy konferentsii "PaVT'2013". Chelyabinsk, 2013, pp. 203-214.

10. Lukin V.V., Marchevskiy I.K., Moreva V.S., Popov A.Yu., Shapo-valov K.L., Shheglov G.A. Uchebno-eksperimentalnyy vychislitel'nyy klas-ter. Primery resheniya zadach [Computing cluster for training and experiments. Part 2. Examples of solving problems]. Vestnik Moskovskogo gosudarstvennogo tekhnicheskogo universiteta imeni N.E. Baumana. Estestvennye nauki, 2012, no. 4, pp. 82-102.

11. Gafurov M.B., Ilgamov M.A. Izgib tsilindricheskoy obolochki konechnoy dliny pri ee poperechnom obtekanii zhidkostyu [Bending of a cylindrical shell of finite length at its crossflow in the fluid]. Prikladnaya mekhanika, 1978, vol. 14, no. 3, pp. 60-67.

12. Ermakov A.V., Shheglov G.A. Modelirovanie metodom vikhrevykh ehlementov dinamiki tsilindricheskoy obolochki v prostranstven-nom potoke zhidkosti [Vortex element method application for a cylindrical shell dynamics simulation in the spatial flow]. Izvestiya vysshikh uchebnykh zavedeniy. Mashinostroenie, 2014, no. 3, pp. 35-41.

Об авторах

Ермаков Андрей Васильевич (Москва, Россия) - аспирант кафедры «Аэрокосмические системы» МГТУ им. Н.Э. Баумана (105005, г. Москва, ул. 2-я Бауманская, 5, e-mail: anvaerm@gmail.com).

Марчевский Илья Константинович (Москва, Россия) - кандидат физико-математических наук, доцент кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана (105005, г. Москва, ул. 2-я Бауманская, 5, e-mail: iliamarchevsky@mail.ru).

Щеглов Георгий Александрович (Москва, Россия) - доктор технических наук, профессор кафедры «Аэрокосмические системы» МГТУ им. Н.Э. Баумана (105005, г. Москва, ул. 2-я Бауманская, 5, e-mail: georg@energomen.ru).

About the authors

Andrey V. Ermakov (Moscow, Russian Federation) - Doctoral Student, Department of Aerospace Systems, Bauman Moscow State Technical University (5, 2nd Baumanskaya st., Moscow, 105005, Russian Federation, e-mail: anvaerm@gmail.com).

Ilya K. Marchevskiy (Moscow, Russian Federation) - Ph. D. in Physics and Mathematical Sciences, Associate Professor, Department of Applied Mathematics, Bauman Moscow State Technical University (5, 2nd Baumanskaya st., Moscow, 105005, Russian Federation, e-mail: iliamarchevsky@mail.ru).

Georgiy A. Shcheglov (Moscow, Russian Federation) - Doctor of Technical Sciences, Professor, Department of Aerospace Systems, Bauman Moscow State Technical University (5, 2nd Baumanskaya st., Moscow, 105005, Russian Federation, e-mail: georg@energomen.ru).

Получено 1.10.2014