Исследование устойчивости конических оболочек вращения с внутренним течением жидкости Текст научной статьи по специальности «Физика»

УДК 539.3+534.1

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

ЖИДКОСТИ

© 2008 С.А. Бочкарев,1 В.П. Матвеенко2

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

Работа поддержана грантом РФФИ (проект №07-08-12144-офи)

Ключевые слова: устойчивость, оболочка вращения, течение жидкости, метод Бубнова-Галеркина, вариационный принцип.

Введение

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

1 Бочкарев Сергей Аркадьевич ([email protected]), Институт механики сплошных сред Уральского отделения РАН, 614013, Россия, г. Пермь, ул. Акад. Королева, 1.

2Матвеенко Валерий Павлович ([email protected]), Институт механики сплошных сред Уральского отделения РАН, 614013, Россия, г. Пермь, ул. Акад. Королева, 1.

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

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

1. Постановка задачи и основные соотношения

Рассматривается упругая коническая оболочка вращения длиной Ь, наименьшим радиусом Я, толщиной Н и углом при вершине а. Внутри оболочки находится идеальная сжимаемая жидкость, которая течет со скоростью и. Требуется найти такую скорость потока, при которой оболочка теряет устойчивость.

Для оболочки приняты гипотезы Кирхгофа-Лява, согласно которым компоненты вектора деформации срединной поверхности, изменения кривизн и кручения в координатной системе записываются следующим образом [8]

du 1 dv u sin a w cos a 1 du dv v sin a

ei = T-. e2 = +-+-' ei2 = + я---'

d s r дв r r r дв d s r

d2w cos a dv 1 d2w sinadw (11)

ds2' r2 дв r d62 r дs'

sin a dv v sin a cos a 1 d2w sin a dw

K12 =--------1---.

r ds r2 r dsde r2 дв

Здесь u, v, w — меридиональная, окружная и нормальная составляющие вектора перемещений, r — текущий радиус оболочки.

Физические соотношения, устанавливающие связь между вектором обобщенных усилий и моментов T и вектором обобщенных деформаций £ = = (ei, Ё2, Ei2,Ki,К2,2ki2}T, представляются в матричном виде

T = (Til, T22, T12, M11, M22, Mi2}T = D£, (1.2)

где ненулевые элементы матрицы D для изотропного материала

D12 = D21 = vDii = VD22 = vEh/(l - v2), D45 = D54 = VD44 = vD55 = vEh2/[l2(l - v2)], D66 = Dssh^12 = Gh3/l2.

Здесь E, v, G — модуль упругости, коэффициент Пуассона и модуль сдвига, соответственно.

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

/вЛТЙ +/6^тр0*5 -/Ь<,тР„5 = 0. (1.3)

5 5 5

Здесь £, Т, З, Р — соответственно векторы обобщенных деформаций, обобщенных усилий и моментов, перемещений, поверхностных нагрузок,

ро = ^ рш<3г,

н

рт —удельная плотность материала оболочки.

Движение идеальной сжимаемой жидкости, расположенной внутри оболочки, и занимающей объем , в случае потенциального течения описывается волновым уравнением, которое в цилиндрических координатах (г, 6, х) имеет вид [9]:

у2 _ <52Ф 1 <92ф <92ф 1 <9ф _ 1 д д 12 дг2 г2 дб2 дх2 г дг е2 дг дх

ф, (1.4)

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

^ = + на Ба = Б/ПБ1. (1.5)

Здесь рf — удельная плотность жидкости; 5 — меридиональная координата оболочки; 5f, 55 — площади, ограничивающие объемы жидкости и оболочки, соответственно. На поверхности раздела оболочка—жидкость 5 а задается условие непроницаемости

дф дw дw

-Г = -г + и-Г на !-6

дп дг д 5

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

ф = 0 при х = 0 и дф/дх = 0 при х = Ь. (1.7)

Применение метода Бубнова-Галеркина [10] к уравнению в частных производных для потенциала возмущения скорости (1.4) с граничными условиями (1.6), (1.7) позволяет получить следующее интегральное соотношение [11]

Шф

1=1

I

dFi dFk 1 dFi dFk / лЛ\ + --r^—rrr + (1 - м )

2\ dFi dFk

"Til — 1V1 ) --

dr dr r2 дб дб ^ ' dx dx

dV

Шф Ш

l=1

Г 2 U dFi

J с2 dx

FkdV

Wf

Шф

ф al +

l=1

fe2

FiFkdV

Фа1 +

ф al-

Wf

(1.8)

z

i=1

¡M-

FkdS

-ai- 2 i=1

f

dNw U—^FkdS д s

Wai = 0, k = 1, Шф.

Здесь Шф, ms — число конечных элементов, на которые разбиваются области жидкости Vf и оболочки Vs; фа, wai — узловые значения для жидкости и оболочки; M = U/c — число Маха; F, N- —функции формы для потенциала возмущений скорости и нормальной составляющей вектора перемещения оболочки.

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

Применяя стандартные процедуры МКЭ к (1.3) и (1.8) получим системы уравнений, которые в объединенном виде можно записать следующим образом

K

d Фа

+ M

d Фа

+ Pf с

d Ф a

+ Pf А

d Фа

= 0,

(1.9)

где

K =

Ks =

Кф =

Сф =

Аф =

Ks

0

0 -P f Кф

, M =

Ms

0

0 -P f Мф

, с =

0 Сф1

ф I a =

Cr^c I '

ф сф

0 As Аф Аф]

шф Vf

2 f F

ms S S o

? V

BTDBdS, Ms =

и

NTP0NdS, As = ^ С UN-dF/дsdS,

í dFT dF 1 dFT dF dFT dF

+

+

dr dr r2 дб дб дх дх

o

dV, Мф = ^ Г4рТрdV>

Шф Vf c

Шф V fV

uF^ds, a; =

д s

Шф

шф Vf

M

2<9FT<9F

д s д s

dV.

Здесь В — матрица связи вектора деформаций £ с вектором узловых перемещений оболочечного конечного элемента; N — матрица функций формы оболочечного элемента.

Традиционно матрица К называется матрицей жесткости, М — матрицей масс, С — матрицей аэродинамического демпфирования и А — матрицей аэродинамической жесткости.

Представляя возмущенное движение оболочки и жидкости следующим образом

Решение поставленной задачи о динамическом поведении оболочек вращения с жидкостью сводится к вычислению и анализу собственных значений X системы (1.10). Для неподвижной жидкости (А = Сф = 0) собственные значения являются действительными. При наличии скорости потока и > 0 собственные значения являются действительными или комплексными. Для оболочки возможны два типа потери устойчивости: статический (дивергенция) и динамический (флаттер). Первый случай неустойчивости характеризуется появлением у одного из собственных значений нулевой действительной части. Во втором случае неустойчивость проявляется в слиянии двух форм колебаний и появлении отрицательной мнимой части у одного из этих собственных значений.

Для вычисления комплексных собственных значений системы (1.10) используется алгоритм на основе метода Мюллера [12].

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

2. Численные примеры

В качестве тестового примера рассматриваются собственные колебания конической оболочки (Е = 6.77 X 1010 н/м2, V = 0.29, рт = 2648 кг/м3, Я = = 0.15 м, Ь = 0.56 м, Н = 5.3 X 10-4 м, а = 15), жестко закрепленной с обоих торцов (и = V = w = дф/д5 = 0) и полностью заполненной неподвижной жидкостью. В табл. 2. представлены низшие собственные частоты колебаний для различных номеров гармоник: в первой строке результаты данной работы; во второй строке результаты конечно-элементных расчетов [13];

З = q ехр(г*Хг), фа = ф ехр(г* Хг),

где q, ф — некоторые функции координат, ? = а X = Х1 + г*Х2 — харак-

теристический показатель, окончательно получим

(1.10)

следующие две строки соответствуют теоретическим и экспериментальным результатам [14]. Полученные результаты хорошо согласуются с известными численными и экспериментальными данными. Расчеты были выполнены при 40 элементах для оболочки и 1000 — для жидкости.

Таблица 1

Собственные частоты колебаний (Гц) конической оболочки, полностью заполненной несжимаемой жидкостью

Номер гармоники j

3 4 5 6 7 8 9

100.86 78.7 63.55 54.23 50.52 52.24 58.20

96.34 75.5 61.07 53.22 50.12 52.14 58.21

101.0 78.7 63.60 54.4 50.80 52.80 —

100.0 76.0 — — 51.00 54.00 —

Далее приведены результаты сопоставления с данными расчетов из работы [7]. Оболочка выполнена из углеродистой стали: средний радиус a = = 0.876 м, L = 0.9144 м, h = 1.5 X 10-3 м, рf = 1000 кг/м3, с = 1500 м/с. Входной и выходной радиусы конических оболочек задаются таким образом, чтобы при любом угле конусности отношения длины оболочки к среднему радиусу и среднего радиуса к толщине оболочки было неизменным, т.е. Lia = 1.0438, a/h = 584. Такой выбор геометрических параметров позволяет сопоставить границы динамической устойчивости конических оболочек с эквивалентной им цилиндрической оболочкой.

В расчетах предполагается, что поток жидкости движется внутри оболочки с постоянной скоростью, так, что U = Q/S [7], где Q — расход жидкости (м3/с) и S —площадь поперечного сечения оболочки.

Рассматривается вариант оболочки свободно опертой (v = w = 0) с двух торцов. На рис. 1а представлены собственные частоты колебаний пустой оболочки (f, рад/с) и оболочки, полностью заполненной жидкостью (ш, Гц), при а = 10 для различных номеров гармоник j. Полученные результаты (для жидкости использовалось 1800 конечных элементов) хорошо согласуются с данными из работы [7]. На рис. 1б для этой же оболочки приводятся значения критической скорости (U, м/с). В этом случае наблюдается количественное расхождение с результатами, полученными в [7]. Более полный набор результатов расчета для рассматриваемого варианта оболочки приведен в табл. 2.

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

Рис. 1. Собственные частоты (а) и критические скорости (б) конической оболочки, свободно опертой с двух торцов, для различных номеров гармоник

Таблица 2

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

] Критическая скорость С/ (м/с)

а = 0 а= = 10 а = 30

[7] Расчет [7] Расчет [7] Расчет

1 147.67 130.23 133.98 122.95 89.59 102.04

3 146.84 129.68 133.15 122.27 89.57 88.81

4 145.6 129.19 131.9 121.59 88.33 87.46

7 87.11 92.09 81.30 88.88 47.29 62.66

9 68.86 71.01 63.87 68.99 38.15 52.61

10 65.54 66.85 61.39 65.39 38.98 53.48

15 82.55 77.11 78.81 76.21 56.41 62.77

Проиллюстрируем обнаруженный феномен на конкретном примере, где для этой же оболочки (] = 1, а = 30) построено изменение действительной и мнимой частей безразмерных собственных значений (П = ~КК (рт/Е)°'5 X 100), соответствующих формам колебаний с числом т — полуволн в меридиональном направлении (рис. 2). Приведенные результаты демонстрируют, что с ростом значения Q первым нулевой величины достигает действительная часть собственного значения, соответствующей форме колебаний с т = 2. Однако при этом соответствующие мнимые части остаются большими нуля и потеря устойчивости не происходит. При дальнейшем увеличении Q обращается в нуль действительная часть собственного значения, соответствующей форме колебаний с одной полуволной. При этом одна из мнимых частей достигает отрицательных значений, т.е. реализуется вариант потери устойчивости, соответствующий дивергенции.

Одновременный анализ картины эволюции действительных и мнимых частей собственных значений позволил установить еще одну особенность.

Reí

12 9 6 3

0 . 0 150 175 200 225 Q 0 150 175 200 225 Q

Рис. 2. Изменения действительной (а) и мнимой (б) частей трех первых безразмерных собственных значений О от расхода жидкости Q для конической оболочки, свободно опертой с двух торцов: линии — расчет; символы—[7]

V. а.

rv* 1 ^ .ж N г - 3 /

m = 2 К:- N ■ N i \ " / / /

1—1-/ \ - \ / >\"

'- к!»

-1 ni - i - [7]

Im

1.2 0.8 0.4 0

-0 4

р-Ф m - 2 л б. « "7

- 3 п к 11 к i ■ и i«

V 'Л ■ г

' " 'А п

m = г 1 >

-

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

Проанализируем устойчивость рассматриваемой оболочки в случае жесткого закрепления торцов. Аналогично предыдущему варианту на рис. 3 представлены собственные частоты колебаний пустой оболочки, оболочки с жидкостью и критические скорости для различных номеров гармоник. В табл. 2. дается более полный набор численных результатов. Как и в предыдущем случае, различие в значениях критических скоростей могут быть объяснены нахождением новых, по отношению к [7], узких зон неустойчивости дивергентного типа. При детальном анализе эволюции действительных и мнимых частей собственных значений (рис. 4, ] = 15, а = 30) также наблюдается существенно немотонное на небольшом интервале изменение мнимых частей (на рис. 4б не показано) при совпадении действительных частей. Обнаруженная закономерность указывает на необходимость крайне аккуратных вычислений, т.к. при не учете возможности такого поведения мнимых частей собственных значений могут быть потеряны некоторые режимы неустойчивости. В частности, для рассматриваемого варианта оболочки при т = 6 (изменение действительной части не представлено) имеет место потеря устойчивости типа флаттер (флаттер по одной степени свободы [15]). Отметим, что аналогичное поведение мнимых частей, обусловленное наличием гидродинамического демпфирования, впервые было

представлено в [16] для случая консольно закрепленной нагруженной цилиндрической оболочки.

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

Таблица 3

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

] Критическая скорость С/ (м/с)

а = 0 а= = 10 а = 30

[7] Расчет [7] Расчет [7] Расчет

1 148.5 130.88 138.96 147.71 96.65 109.5

3 147.26 130.31 137.3 145.55 96.65 116.13

5 145.6 129.12 134.4 144.54 91.25 117.75

7 124.03 126.29 116.15 121.84 68.85 90.18

10 96.65 95.24 89.59 92.11 53.92 68.40

12 91.26 87.41 85.45 85.10 54.75 80.16

15 92.92 85.66 88.77 84.26 63.45 85.21

Результаты численных экспериментов по исследованию устойчивости конических оболочек взаимодействующих с потоком жидкости приводят и к другим качественным выводам. Первый из них заключается в том, что при увеличении угла конусности имеет место сгущение спектра собственных частот колебаний оболочки, как с неподвижной, так и текущей жидкостью. Иллюстрацией к данному выводу является приведенные на рис. 5а расчеты (1000 элементов для жидкости) для оболочки со следующими параметрами: Е = 2.06 X 1011 н/м2, V = 0.3, рш = 7850 кг/м3, Я = 0.2 м, Ь = 0.4 м, Н = 2X10-3 м (радиус входного сечения остается неизменным, а при изменении угла конусности меняется длина образующей оболочки). Этот вывод с учетом предыдущих результатов имеет важное значение для выбора страте-

Рис. 4. Изменения действительной (а) и мнимой (б) частей безразмерных собственных значений О от расхода жидкости Q для конической оболочки, жестко закрепленной на обоих торцах: линии — расчет; символы — [7]

гии вычислений, связанный с обеспечением детального анализа поведения действительных и мнимых частей собственных значений.

Рис. 5. Собственные частоты колебаний ю (Гц) оболочки с жидкостью (а) и критические скорости и (м/с) жидкости (б) при различных углах конусности а

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

Заключение

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

Литература

[1] Lakis, A.A. Dynamic analysis of anisotropic fluid-filled conical shells / A.A.Lakis, P. Van Dyke, H. Ouriche // J. Fluids Struct. - 1992. - V. 6. -№2. - P. 135-162.

[2] Selmane, A. Vibration analysis of anisotropic open cylindrical shells subjected to a flowing fluid / A. Selmane, A.A. Lakis //J. Fluids Struct. -1997. - V. 11. - P. 111-134.

[3] Zhang, Y.L. A finite element method for modelling the vibration of initially tensioned thin-walled orthotropic cylindrical tubes conveying fluid / Y.L.Zhang, D.G.Gorman, J.M.Reese // J. Sound Vib. - 2001. - V. 245. -№1. - P. 93-112.

[4] Zhang, Y.L. Vibration of prestressed thin cylindrical shells conveying fluid / Y.L.Zhang, D.G.Gorman, J.M.Reese // Thin-Walled Struct. -2003. - V. 41. - P. 1103-1127.

[5] Kochupillai, J. A semi-analytical coupled finite element formulation for shells conveying fluids / J. Kochupillai, N. Ganesan, C. Padmanabhan // Comput. Struct. - 2002. - V. 80. - P. 271-286.

[6] Kochupillai, J. A semi-analytical coupled finite element formulation for composite shells conveying fluids / J. Kochupillai, N. Ganesan,

C. Padmanabhan // J. Sound Vib. - 2002. - V. 258. - №2. - P. 287-307.

[7] Kumar, D.S. Dynamic analysis of conical shells conveying fluid /

D.S.Kumar, N. Ganesan // J. Sound Vib. - 2008. - V. 310. - №1-2. -P. 38-57.

[8] Бидерман, В.Л. Механика тонкостенных конструкций / В.Л. Бидерман. - М.: Машиностроение, 1977. - 488 с.

[9] Вольмир, А.С. Оболочки в потоке жидкости и газа. Задачи гидроупругости / А.С. Вольмир. - М.: Наука, 1979. 320 с.

[10] Флетчер, К. Численные методы на основе метода Галеркина / К. Флет-чер. - М.: Мир, 1988. - 352 с.

[11] Бочкарев, С.А. Численное исследование влияния граничных условий на динамику поведения цилиндрической оболочки с протекающей жидкостью / С.А. Бочкарев, В.П. Матвеенко // Известия РАН. МТТ. -2008. - №3. - С. 190-201.

[12] Матвеенко, В.П. Об одном алгоритме решения задачи о собственных колебаниях упругих тел методом конечных элементов / В.П. Матвеенко // Краевые задачи теории упругости и вязкоупругости. - Свердловск, 1980. - С. 20-24.

[13] Григорьев, В.Г. Методология исследования динамических свойств сложных упругих и гидроупругих систем: дис. ... д-ра техн. наук. /

B.Г. Григорьев. - М., 2000.

[14] Горбунов, Ю.А. Теоретическое и экспериментальное исследование спектра собственных неосесимметричных колебаний конической оболочки с жидкостью при наличии внутреннего давления / Ю.А.Горбунов, Л.М. Новохатская, В.П.Шмаков // Динамика упругих и твердых тел, взаимодействующих с жидкостью. - Томск: Изд-во ТГУ, 1975. -

C. 47-52.

[15] Paidoussis, M.P. Flutter of thin cylindrical shells conveying fluid / M.P. Paidoussis, J.-P. Denise // J. Sound Vib. - 1972. - V. 20. - №1. -P. 9-26.

[16] Горачек, Я. Влияние закрепления краев цилиндрической оболочки с протекающей жидкостью на ее динамические характеристики / Я. Горачек, И. Золотарев // Прикл. механика. - 1984. - Т. 20. - №8. -С. 88-98.

Поступила в редакцию 15/VII/2008;

в окончательном варианте — 15/ VII/2008.

AN INVESTIGATION OF INTERNAL STABILITY OF CONICAL SHELLS WITH INTERNAL FLUID FLOW

© 2008 S.A. Bochkarev, V.P. Matveyenko3

A mixed finite-element algorithm has been proposed for studying dynamic stability of the conical shells of revolution interacting with the internal flow of the compressible fluid. The fluid motion is described in the framework of the potential theory, whose equations are reduced to the integral form by means of the Bubnov-Galerkin method. Treatment of the shell is based on the variational principle of virtual displacements, which includes a linearized Bernoulli equation for calculation of hydrody-namic pressure acting on the shell on the side of the fluid. Solution of the problem reduces to evaluation and analysis of the eigenvalues of the connected systems of equations. As a numerical example, the problem on the influence of the cone angle on the boundary of the shell hydroelastic stability is shown.

Keywords and phrases: stability, shell of revolution, liquid flow, Bubnov-Galerkin method, variational principle.

Paper received 15/V77/2008. Paper accepted 15/V77/2008.

3 Bochkarev Sergey Arkadeivich ([email protected]), Matveyenko Valery Pavlovich ([email protected]), Institute of Continuous Media Mechanics of the Ural Branch of the Russian Academy of Sciences, Perm, 614013, Russia.