МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ТЕЧЕНИЯ ВЯЗКО-ПЛАСТИЧЕСКОЙ ЖИДКОСТИ В СМЕШАННЫХ ПЕРЕМЕННЫХ "ВИХРЬ-СКОРОСТЬ" Текст научной статьи по специальности «Физика»

УДК 519.63

DOI: 10.18698/2309-3684-2021-2-314

Математическое моделирование течения вязко-пластической жидкости в смешанных переменных «вихрь-скорость»

© С.В. Лебедев, С.Р. Абдуллин, Н.И. Бондаренко МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

В работе рассматривается задача о численном трехмерном моделировании движения вязко-пластических жидкостей Шведова-Бингама в цилиндрических установках. Такие задачи имеют практическое значения для ряда отраслей промышленности: химической, пищевой и некоторых, при создании перспективных установок с перемешиванием и/или обработкой жидких сред. Жидкости в данных установках проявляют, как правило, сильно вязкие свойства, которые наиболее адекватно описываются моделями неньютоновских вязкопластических сред. Неньютоновские свойства жидкостей усложняют уравнения движения, что приводит к необходимости разработки метода численного расчета данных уравнений стационарного пространственного течения жидкой среды. В данной работе предложен численный алгоритм решения трехмерной задачи движения жидкости Шведова-Бингама, основанный на использовании переменных «вихрь-скорость» и применении неявной конечно-разностной схемы расчета. Представлены некоторые результаты расчета трехмерного движения жидкости Шведова-Бингама в цилиндрической трубе.

Ключевые слова: неньютоновские жидкости, вязко-пластическая среда Шведова-Бингама, нелинейная вязкость, неявная конечно-разностная схема, переменные «вихрь-скорость»

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

Различные модели нелинейных реологических уравнений неньютоновских жидкостей предложены в работах [1-13]. Возможен учет предыстории напряжения и деформации, например, на основе моделей Максвелла, Фойгта и Кельвина [1].

В данной работе рассматривается модель вязко-пластической жидкости Шведова-Бингама [2, 4-10], которая в силу специфики

сильно вязких жидких сред достаточно адекватно описывает движение жидкостей в рассматриваемых перспективных установках. Жидкая среда Шведова-Бингама имеет начальный предел текучести, ниже которого она не течет и имеет свойства твердого тела.

Реологическая модель вязко-пластической жидкой среды. В трехмерной (пространственной) тензорно-инвариантной постановке вязко-пластическая модель Шведова-Бингама выражается следующими формулами:

Б - 0, т{2)<т0, (1а)

т = (2^0/^22)т(2)~то, (1б)

где Б — тензор скоростей деформации, т — тензор-девиатор напряжений, компоненты которого выражаются по формулам т - а -1/3 рЕ ; о — тензор напряжений, р — давление в жидкости, ^(2), т(2) — вторые инварианты (интенсивности) [14] тензоров Б и т, т0 — предельное напряжение сдвига (начальный предел текучести) [4], л — коэффициент вязкости.

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

т = 2//тахБ, (2)

где //^//^»1 (обычно достаточно Лтж/Мо ~ 12-15 ) (при большой

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

меняется от точки к точке в соответствии с установленными выше соотношениями (1)-(2).

Обобщая (1), (2), имеем реологическое уравнение для нелинейно-вязкой модели течения:

т - Б.

Система уравнений течения несжимаемой вязко-пластической жидкой среды. Математическая постановка задачи движения

рассматриваемой несжимаемой жидкой среды включает уравнение неразрывности и уравнение движения [11]:

V- и = 0, (3)

0 1

— + и-V® и---Vp + V■ т, (4)

дг р

где и — вектор скорости, р — плотность, р — давление жидкой среды, V —набла-оператор [14], • —знак скалярного произведения векторов, ® — знак тензорного умножения [14].

К рассматриваемой математической модели следует присоединить уравнения (1)-(2).

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

Вихрь (завихренность) скорости определяется как ротор вектора скорости [3,11]: % = Vxи , здесь х — знак векторного произведения. Вихрь характеризует степень закрученности потока в каждой точке. Чтобы получить гидродинамические уравнения (3), (4) в смешанных переменных, применим к обеим частям уравнения (4) дифференциальный оператор ротор.

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

+ Vх(% х и + % ^у^ ), (5)

% ^х и, (6)

V-и = 0, (7)

где тензор скоростей деформаций выражается по следующим формулам:

Б = + V® уг). (8)

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

Алгоритм численного решения системы уравнений течения вязко-пластической жидкости Для расчетов стационарного течения в соответствии с уравнениями (5)-(8) используется метод установления, таким образом, расчет производится до «обнуления» частной производной по времени в (5).

Расчет вектора вихря и вектора скорости на последующем временном шаге проводится в два шага:

• по (5) рассчитываются компоненты вектора вихря;

• по (6), (7) рассчитываются компоненты вектора скорости (с учетом изменения вектора вихря).

Алгоритм расчета компонент вектора вихря (Д)"+1 • Рассмотрим ^ (I = 1,2,3) — криволинейные ортогональные координаты, и обозначим далее щ, Д — компоненты векторов скорости и вихря в этой криволинейной ортогональной системе координат, ^ — коэффициенты Ляме данной криволинейной системы координат.

Для расчета трехмерного вектора вихря на последующем временном шаге в соответствии с уравнением (5) применяется разработанная двухэтапная конечно-разностная схема второго порядка точности аппроксимации по пространственному и временному шагам.

На первом этапе расчета изменений трех компонент вектора вихря

*3* / \*3* / \*3*

, (ОД) и (ОД) используется неявная многошаговая схема

Освальда, основанная на одной из модификаций метода переменных направлений и состоящая из девяти подэтапов.

1) «Явный» этап для компоненты вихря Д :

(ОД Г = А/ . {[ В, ] ( Д )" + [ ^3 1 ( Д )" - [ + А ] (Д )" - А } -

2) Первый «неявный» этап для компоненты вихря Д (с учетом изменения Д ):

(1+ А// 2.)(ОД Г = = Аt ■ {[В21 ] ( Д )" + [В23 ] ( Д )" - [Л21 + А23 ] ( Д)" - Ц} +

+ А//2■[В21 ](ОДГ.

3) Первый «неявный» этап для компоненты вихря Д (с учетом изменения Д ):

(1+ А/ 2 ■ А31 )(ОД Г =

= А/ ■ {[ В31 ] (Д )" + [ В32 ] (Д)" - [ 4, + А32 ] (Д)" - А} +

+ А//2 ■[В31 ](ОД )*".

4) (ОД Г=(Д )'1' -

5) Первый «неявный» этап для компоненты вихря Д (с учетом изменения Д ):

(1 + Лг/ 2 ■ А2 )8 Г= (йГ+Лг/2 ■[ 512 ](й )".

6) Второй «неявный» этап для компоненты вихря й (с учетом изменения ):

(1 + Лг/ 2 ■ А32) (83 Г = 8 Г + Лг/ 2 ■ [ ^32 ] 8 Г.

7) (83 ) = 8з Г.

8) Второй «неявный» этап для компоненты вихря (с учетом изменения ):

(1 + Лг/ 2 ■ 4з) 8 Г = 8 Г + лг/ 2 ■ [ Д23 ] 8 Г.

В используемой схеме конечно-разностные операторы А и в

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

1

д

И2 Л3 [д^2 — -1—

—-1-д-—

КК 1^3

К3 -уеЖ

%3иА

К д(£К)

КК2 д^2 ,

К2 д(йЛ)"

КК д^3 _

К д(Й2К )

¿А 1

К2 д(£К)

о А,

о А

13'

о Д

Д13.

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

8 Г =Лг-{[ Д2 ]((й)" + (82 )*372) + [ Д13 ]((й)" + 8 Г/2

-Лг-{[ А + А13 ]((?)" +(й )*372) +А},

8 Г = лг - {[Д23 ] ((й)" + 8 Г72)+[В21 ] ((й)" + (8?, Г72

{[ А23 + А,! ]((й)" + (82 ГД) + В2},

-Лг -<

(ОД )"+1 =А/ ■{[ В3, ]((Д)" + (ОД Г/2 ) + + [ В32 ]((Д )" + (ОД Г/2)}-

■{[ 4, + А, ]((Д)" +(Д ГД)+А}.

-А/ ■{[ 4, + А31 ]((Д)" +

Алгоритм расчета компонент вектора скорости (щ )"+1. Задача

определения поля скорости по тождественно равной нулю дивергенции и по известному полю ротора (уравнения (6), (7)) получило название задачи Коши-Римана. Для численного решения задачи Коши-Римана в конечных разностях важно обеспечить равенство количества переменных (неизвестных значений скорости) и количества независимых линейных уравнений. Данное условие обеспечивается при использовании определенной разнесенной трехмерной конечно-разностной сетки, в которой компоненты вихря и скорости определяются в несовпадающих точках.

Рассмотрим элементарный расчетный объем в криволинейной ортогональной системе координат дг (г = 1,2,3), имеющий форму

прямоугольного параллелепипеда; определяем каждую из трех компонент скорости в центре грани, перпендикулярной компоненте скорости, каждую из трех компонент вектора вихря скорости — в середине ребра, параллельной компоненте вихря. Благодаря данному расположению:

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

• конечно-разностные аналоги уравнений (6) и (7) просто и наглядно аппроксимируют уравнения со вторым порядком точности по пространственным переменным — конечно-разностные шаблоны уравнений представлены на рис. 1 (определение вихря скорости) и на рис. 2 (уравнение неразрывности).

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

Рис. 1. Шаблон определения компонент вихря (синие стрелки — компоненты скорости, красная стрелка — компонента вихря)

Рис. 2. Шаблон для уравнения неразрывности (синие стрелки — компоненты скорости)

Рис. 3. Иллюстрация метода ортогональных проекций (на примере двух уравнений с двумя неизвестными)

Граничные и начальные условия. Граничные условия для скорости: на входе в рабочую область задается постоянный профиль

и

2

и

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

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

$*= (Ух«I ,

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

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

и

2+1/2

(з< - <-1 V2, (г11 - 2 (С )„-(¿Г1 I, (* -1,2,3).

Метод апробирован на двух примерах (получено удовлетворительное совпадение результатов):

• течение в кубе с подвижной верхней крышкой;

• течение в слое между двумя концентрическими цилиндрами. Результаты численного моделирования течения в цилиндрической трубе. Проведена серия расчетов жидкости Шведова-Бингама в цилиндрической трубе; в данной серии варьируется входная скорость (по сравнению с «базовой» исходной скоростью (и0 - 2 см/с ).

Все расчеты сравниваются с данным «базовым» расчетом, в котором радиальный размер ядра течения составляет 60% от радиуса трубы. Характерный профиль скорости представлен на рис. 4.

г 1,00

0,75

0,50 0,25 0,00 0,25 0,50 0,75 1,00

0 12 3 4

0

1

2

3

4

Рис. 4. Профиль осевой скорости вязко-пластического течения в цилиндрической трубе

г

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

Таблица 1

Влияние скорости (на входе) вязко-пластического течения в трубе на размеры твердого ядра (в относительных величинах)

0,5 1,0 1,5 2,0 3,0 4,0 5,0

в/R 0,70 0,60 0,51 0,45 0,38 0,32 0,28

Выводы. Разработана конечно-разностная схема для моделирования пространственных течений вязких и неньютоновских жидких сред. На основе анализа реологических характеристик жидких сред в рассматриваемых перспективных устройствах перемещения и обработки для расчетов предложена реологическая модель среды с начальным пределом текучести — вязко-пластическая жидкость Шведова-Бингама. Для численных расчетов разработана численная модификация реологической модели с введением максимально возможной вязкости. Проведены расчеты вязко-пластического течения в круглой трубе, найдены зависимости размеров ядра течения от входной скорости, что определяет диапазон возможностей при выборе оптимальных режимов работы установок. Разработанный метод расчета будет востребован и для расчета течений в аналогичных установках со сложной геометрией.

ЛИТЕРАТУРА

[1] Астарита Дж. Основы гидромеханики неньютоновских жидкостей. Москва, Мир, 1978, 312 с.

[2] Климов Д.М., Петров А.Г., Георгиевский Д.В. Вязкопластические течения: динамический хаос, устойчивость, перемешивание. Москва, Наука, 2005, 394 с.

[3] Димитриенко Ю.И. Нелинейная механика сплошной среды. Москва, Физма-тлит, 2009, 624 с.

[4] Матвиенко О.В., Базуев В.П., Сабылина Н.Р., Асеева А.Е., Суртаева А.А. Исследование установившегося течения вязкопластического битумного вяжущего, описываемого моделью Шведова-Бингама, в цилиндрической трубе. Вестник Томского государственного архитектурно-строительного университета, 2019, т. 21, № 3, с. 158-177.

[5] Кийко И.А. Вязко-пластическое течение материалов. Физико-математические основы технологии обработки давлением. Москва, Механико-математический факультет МГУ, 2001, 136 с.

[6] Евтюков С.А., Овчаров А.А., Замараев И.В. Построение механореологиче-ских моделей процессов взаимодействия рабочих органов строительно-дорожных машин со средой: учебное пособие. Санкт-Петербург, СПбГАСУ, 2011, 59 c.

[7] Трусов П.В., Швейкин А.И. Теория пластичности: учебное пособие. Пермь, Изд-во Пермского нац. исследовательского политехнического ун-та, 2011, 418 c.

[8] Cheddadi I., Saramito P., Dollet B., Raufaste C., Graner F. Understanding and predicting viscous, elastic, plastic flows. European Physical Journal E, 2011, vol. 34, iss. 1, art. no. 1. DOI: 10.1140/epje/i2011-11001-4

[9] Oleinik E.F. Plasticity of semicrystalline flexible-chain polymers at the microscopic and mesoscopic levels. Polymer Science — Series C, 2003, vol. 45, iss. 1, pp. 17-117.

[10] Drozdov A.D., Klitkou R., Christiansen J.D. Cyclic viscoplasticity of semicrystalline polymers with finite deformations. Mechanics of Materials, 2013, vol. 56, pp. 53-64.

[11] Кочин Н.Е., Кибель И.А., Розе Н.В. Теоретическая гидромеханика. Ч.2. Москва, Физматгиз, 1963, 728 с.

[12] Dimitrienko Y.I., Li S. Modeling of non-Newtonian resin flows in composite microstructures. IOP Conference Series: Materials Science and Engineering, 2019, vol. 683, iss. 1, art. no. 012008. DOI: 10.1088/1757-899X/683/1/012008

[13] Dimitrienko Y.I., Li S. Modeling the nonlinear permeability of porous composite structures with non-Newtonian fluid. IOP Conference Series: Materials Science and Engineering, 2020, vol. 934, iss. 1, art. no. 012014. DOI: 10.1088/1757-899X/934/1/012014

[14] Dimitrienko Yu.I. Mekhanika sploshnoj sredy. T. 1. Tenzornyj analiz [Continuum Mechanics. Vol. 1. Tensor analysis]. Moscow, BMSTU Publ., 2011, 367 p.

Статья поступила в редакцию 19.04.2021

Ссылку на эту статью просим оформлять следующим образом:

Лебедев С.В., Абдуллин С.Р., Бондаренко Н.И. Математическое моделирование течения вязко-пластической жидкости в смешанных переменных «вихрь-скорость».

Математическое моделирование и численные методы, 2021, № 2, с. 3-14.

Лебедев Сергей Витальевич — канд. физ.-мат. наук, доцент кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана, инженер Научно-образовательного центра «Суперкомпьютерное инженерное моделирование и разработка программных комплексов» МГТУ им. Н.Э. Баумана. e-mail: serleb@rambler.ru

Абдуллин Салават Роальдович — старший преподаватель кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. e-mail: mai-sal@yandex.ru

Бондаренко Николай Иванович — канд. техн. наук, доцент кафедры «Теоретическая механика» МГТУ им. Н.Э. Баумана. e-mail: colia.bond@yandex.ru

Mathematical simulation of viscous-plastic fluid at mixed variables «curl-velocity»

© S.V. Lebedev, S R. Abdullin, N.I. Bondarenko

Bauman Moscow State Technical University, Moscow, 105005, Russia

The paper considers the problem of numerical three-dimensional modeling of the movement of Shvedov-Bingham viscous-plastic fluids in cylindrical installations. Such tasks are of practical importance for a number of industries: chemical, food and some, when creating promising installations with mixing and / or processing of liquid media. Fluids in these installations, as a rule, exhibit highly viscous properties, which are most adequately described by models of non-Newtonian viscoplastic media. The non-Newtonian properties of liquids complicate the equations of motion, which leads to the need to develop a methodfor the numerical calculation of these equations for the stationary spatial flow of a liquid medium. In this paper, we propose a numerical algorithm for solving the three-dimensional Shvedov-Bingham fluid motion problem, based on the use of vortex-velocity variables and the use of an implicit finite-difference calculation scheme. Some results of calculating the three-dimensional motion of the Shvedov-Bingham fluid in a cylindrical tube are presented.

Keywords: non-Newtonian fluids, Shvedov-Bingham visco-plastic medium, implicitfinite-difference scheme, variables «curl-velocity»

REFERENCES

[1] Astarita J. Osnovy gidromekhaniki nen'yutonovskih zhidkostej [Fundamentals of hydromechanics of non-Newtonian fluids]. Moscow, Mir Publ., 1978, 312 p.

[2] Klimov D.M., Petrov A.G., Georgievsky D.V. Vyazkoplasticheskie techeniya: dinamicheskij haos, ustojchivost', peremeshivanie [Viscoplastic flows: dynamic chaos, stability, mixing]. Moscow, Nauka Publ., 2005, 394 p.

[3] Dimitrienko Yu.I. Nelineynaya mekhanika sploshnoy sredy [Nonlinear continuum mechanics]. Moscow, Fizmatlit Publ., 2009, 624 p.

[4] Matvienko O.V., Bazuev V.P., Sabylina N.R., Aseeva A.E., Surtaeva A.A. Shvedov-Bingham model of steady flow of visco-plastic bitumen binder in cylindrical tube. Vestnik of Tomsk State University of Architecture and Building. English version appendix, 2019, vol. 21, no. 3, pp. 158-177.

[5] Kiiko I.A. Vyazko-plasticheskoe techenie materialov. Fiziko-matematicheskie osnovy tekhnologii obrabotki davleniem [Visco-plastic flow of materials. Physical and mathematical foundations of pressure treatment technology]. Moscow, Moscow State University Faculty of Mechanics and Mathematics, 2001, 136 p.

[6] Evtyukov S.A., Ovcharov A.A., Zamaraev I.V. Postroenie mekhanoreologiche-skih modelej processov vzaimodejstviya rabochih organov stroitel'no-dorozhnyh mashin so sredoj: uchebnoeposobie [Construction of mechanorheological models of processes of interaction of working bodies of road construction machines with the environment: textbook]. St. Petersburg, SPbGASU Publ., 2011, 59 p.

[7] Trusov P.V., Shveykin A.I. Teoriya plastichnosti: uchebnoe posobie [Theory of plasticity: textbook]. Perm, Perm National University Research Polytechnic University Publ., 2011, 418 p.

[8] Cheddadi I., Saramito P., Dollet B., Raufaste C., Graner F. Understanding and predicting viscous, elastic, plastic flows. European Physical Journal E, 2011, vol. 34, iss. 1, art. no. 1. DOI: 10.1140/epje/i2011-11001-4

[9] Oleinik E.F. Plasticity of semicrystalline flexible-chain polymers at the microscopic and mesoscopic levels. Polymer Science — Series C, 2003, vol. 45, iss. 1, pp. 17-117.

[10] Drozdov A.D., Klitkou R., Christiansen J.D. Cyclic viscoplasticity of semicrystalline polymers with finite deformations. Mechanics of Materials, 2013, vol. 56, pp. 53-64.

[11] Kochin N.E., Kibel I.A., Rose N.V. Teoreticheskaya gidromekhanika. Ch. 2. [Theoretical Hydromechanics. Part. 2]. Moscow, Fizmatgiz Publ., 1963, 728 p.

[12] Dimitrienko Y.I., Li S. Modeling of non-Newtonian resin flows in composite microstructures. IOP Conference Series: Materials Science and Engineering, 2019, vol. 683, iss. 1, art. no. 012008. DOI: 10.1088/1757-899X/683/1/012008

[13] Dimitrienko Y.I., Li S. Modeling the nonlinear permeability of porous composite structures with non-Newtonian fluid. IOP Conference Series: Materials Science and Engineering, 2020, vol. 934, iss. 1, art. no. 012014. DOI: 10.1088/1757-899X/934/1/012014

[14] Dimitrienko Yu.I. Mekhanika sploshnoj sredy. T. 1. Tenzornyj analiz [Continuum Mechanics. Vol. 1. Tensor analysis]. Moscow, BMSTU Publ., 2011, 367 p.

Lebedec S.V., Cand. Sc. (Phys.—Math.), Assoc. Professor of Department of Computational Mathematics and Mathematical Physics, Bauman Moscow State Technical University, engineer of Research and Education Center Scientific and Educational Centre of Supercomputer Engineering Modeling and Software Packages (Simplex), Bauman Moscow State Technical University. e-mail: serleb@rambler.ru

Abdullin S.R., Assistant Prof. of Department of Computational Mathematics and Mathematical Physics, Bauman Moscow State Technical University. e-mail: mai-sal@yandex.ru

Bondarenko N.I., Cand. Sc. (Eng.), Assoc. Professor of Theoretical Mechanic, Bauman Moscow State Technical University. e-mail: colia.bond@yandex.ru