CC BY
52
УДК 532.517.2:534.2
ИССЛЕДОВАНИЕ НЕНЬЮТОНОВСКОГО ЭФФЕКТА ПРИ ОБТЕКАНИИ ЦИЛИНДРА ПОЛЗУЩИМ ПОТОКОМ УПРУГОВЯЗКОЙ ЖИДКОСТИ В ПЛОСКОМ КАНАЛЕ
© 2009 Б.А. Снигерев1
Рассмотрена модель течения полимерной жидкости в канале вискозиметра с падающим грузом. В качестве реологического конститутивного соотношения использованы модель упруговязкой жидкости Олдройда-Б; релаксационная модель упруговязкой жидкости Фан-Тьен — Таннера (РТТ), предсказывающая аномалию вязкости раствора полимера и наличие продольной вязкости в течениях, обусловленных нормальным напряжением. Для этих моделей получены численные решения задачи о падении груза при различных граничных условиях на стенках канала. Результаты расчетов удовлетворительно согласуются с численными и экспериментальными данными других авторов.
Ключевые слова: вязкоупругая жидкость, вискозиметр, метод конечных элементов.
Введение
Движение тел сферической или цилиндрической формы в жидкостях, обладающих неньютоновскими реологическими свойствами, широко применяется для исследования процесса седиментации при производстве наполненных полимерных систем, а также в вискозиметрии при измерениях вязкости жидкостей, в том числе растворов и расплавов полимеров [1]. Известно, что растворы и расплавы полимеров демонстрируют градиентную зависимость продольной и сдвиговой вязкости, первую и вторую разности нормальных напряжений и эффекты запаздывания, то есть являются нелинейными вязкоупругими средами. Поэтому возникает необходимость более подробного исследования поведения жидкостей с особыми свойствами в каналах вискозиметра с падающим грузом с помощью численных методов,
хСнигерев Борис Александрович (snigerev@mail.knc.ru), лаборатория моделирования технологических процессов института механики и машиностроения КазНЦ РАН, 420111, Россия, г. Казань, ул. Лобачевского, 2/31.
дающих возможность детально описывать сложные неоднородные течения полимерных сред в различных участках вискозиметра.
Согласно общей схеме течения в вискозиметре с падающим грузом, на груз сферической формы действуют следующие силы: сила Архимеда Farh; вес груза и суммарная сила вязкого трения Fvisc. В случае движения груза сферической формы в вязкой жидкости, находящейся в неограниченном контейнере, суммарная сила вязкого трения Fvisc определяется в соответствии с формулой Стокса: Fvisc = 6n^rv. При достижении грузом состояния равномерного и прямолинейного движения v = const эти силы уравновешиваются.
Пусть груз сферической формы имеет радиус r и плотность р и движется под действием сил гравитации, трения и силы Архимеда в цилиндрическом контейнере радиуса R. Упруговязкая жидкость, находящаяся в контейнере, характеризуется плотностью g>, вязкостью при нулевой скорости сдвига По и характеристическим временем релаксации напряжения. Полимерная жидкость может также содержать растворитель c вязкостью ns и обладать эффектом аномалии вязкости.
Условие равновесия груза приводит к следующему уравнению:
4 о , 4 о
Farh + Fvisc + (-mg) = 3nr gg + 6n^rv - 3nr gg = 0.
Откуда скорость падения груза определяется по формуле
2 r2g ( ,)
V = 9Т(g - g )-
Из последней формулы путем измерения скорости v можно оценивать значение вязкости исследуемой жидкости.
Несмотря на кажущуюся простоту, течение жидкости вблизи груза достаточно сложное, а именно течение вблизи передней и задней кромки груза происходит в основном под действием нормальных напряжений, а течение в канале между грузом и стенкой контейнера происходит под действием значительных сдвиговых напряжений. Относительная важность вклада различных факторов, влияющих на скорость падения груза, зависит в первую очередь от реологических свойств исследуемой жидкости, веса груза и относительных размеров груза и контейнера. В настоящей работе получено численное решение двумерной задачи об обтекании бесконечного цилиндра в плоском канале, а текже исследованы две схемы обтекания цилиндра.
1. Постановка задачи
Рассмотрим постановку задачи течения вязкоупругой жидкости в канале при обтекании цилиндра под действием перепада давления. Уравнения движения запишем в виде
/ с)ю \
+ р = -УР + У-т> (1Л)
V- V = 0, (1.2)
где д — плотность жидкости, V — вектор скорости, Р — давление, т — девиатор напряжения.
Реологическое уравнение релаксационного типа запишем в виде
т = т 1 + т 2, (1.3)
где т 1 — неньютоновская составляющая напряжения, определяемая из конститутивного соотношения [3, 5]
т 1 + Л
£ + /1(т иО)
+ ¡2(Л,т 1)т 1 = 2 щП. (1.4)
Ньютоновская составляющая общего напряжения определяется по формует 1
ле т2 = 2^20. В (1.4) —--верхняя конвективная производная от тензо-
Ы
ра напряжения т 1, П1 — вязкость полимера, П2 — вязкость растворителя, Л — время релаксации напряжений. Выбор параметров /1,/2 определяет вид конститутивного реологического соотношения. При /1 =0 и /2 = 0 получается контравариантная модель упруговязкой жидкости Олдройд-Б [4, 5], а выбор функций /1, /2 в виде
/1(тП)= £(Б • т + т • Б), /2(Л,т) = вхр^Г^ 1т(т) (1.5)
определяет экспоненциальную модель Фан-Тьен — Таннера (РТТ) [3, 4]. В уравнении (1.5) £ртт и £ являются материальными параметрами, определяющими поведение данной модели, Ьт(т) = тхх +туу — след тензора. Верхняя конвективная производная определяется как
5т дт чТ
— = — + V • Ут -V— • т — т • (V—)1,
бг дг
тензор скоростей деформации определяется соотношением П = ^ (Vv + + Vvт).
Уравнения (1.1)—(1.4) запишем в безразмерном виде введением следующих переменных:
х* = - у* = - и* = — —* = — Г = — (1б)
х = ь' у = ь' и = и' и = и' г = ь' (1)
ьР ьт
р =
(П1 + П2) и' (П1 + П2)и'
где и— характеристическая скорость задачи, а ь— характерный размер, соответствующий высоте канала. Применение характерных масштабов (1.6) к уравнениям (1.1)—(1.4) приводит к безразмерному виду
д— , „ \ „ „ .
Яе[ — + и •Уи] = —УР + вУ2— + V- т 1, (1.7)
V- и = 0, (1.8)
Ше
т? + /1(т "П)
+ /2 (Л,т 1)т 1 =2(1 — в)Б, (1.9)
где безразмерные комплексы Ее , Ше и в являются числами Рейнольдса, Вайссенберга и коэффициентом ретардации соответственно
Ее = Ше = в =
пг + П2 I* т + П2
2. Метод решения
Для решений уравнений (1.7)—(1.9) применяется метод конечных элементов. Для этого запишем интегральные тождества по методу Бубнова-Галер-кина [6]
/ Ее(дЫ + У ^ = I + + V тг) wdQ,
^ ад dQ = ! ((V • V) д dQ = 0,
/п
Ше
I + «т
+ ¡2(т 1,Б)г^ wdО = ! (1 - в) (Vv + VvT) wdn.
Далее область течения О покрывается семейством четырехугольных элементов Ое. Для скоростей применяется квадратичная аппроксимация ф = = {фг,1 = 1, 9} и линейная аппроксимация ф = {ф^,г = 1, 4} для давления и напряжений. Интегралы на четырехугольнике вычисляются с помощью 4-точечной квадратурной формулы Гаусса. После применения стандартной процедуры метода конечных элементов получаем матричные уравнения следующего вида:
ReMvn+1 + вAvn+1 - Ерп+г = КеМюп - Ее^п + Етп,
= 0,
где М — матрица масс, Е — дискретизация оператора дивергенции, Ет — транспонированная матрица оператора дивергенции, А — сеточный аналог оператора Лапласа, С — дискретизация конвективного члена. Уравнения для напряжении считаются раздельно при известном распределении поля скоростей в виде уравнений следующего вида:
Мтп+1 = Мтп - 0(ип)тп + Е(тп)тп+\
где М — матрица масс, О — дискретизация оператора верхней конвективной производной, Е — дискретизация дополнительных членов уравнения (1.4). Для стабилизации счета применяется дополнительная вязкость вдоль линий тока в виде
N * = N + (XV • VN. (2.1)
Для вычисления интегралов матриц конвективных членов вводится весовая функция N*, определяемая по формуле (2.1). Формулы для вычисления параметра а, зависящего от размера конечного элемента и средней скорости в данном элементе, приведены в работе [9]. Для решения уравнений на
п
каждом временном слое применяется прямой метод разложения Холецко-го [6]. Решение начально-краевой задачи получено с помощью процедуры установления по времени. Для обеспечения сходимости решения по сетке задача решалась на последовательности сгущающихся сеток с числом узлов конечно-элементной сетки N = 12000, 24564. Сетка конечных элементов существенно сгущается вблизи цилиндра для более точного учета концентрации больших напряжений.
3. Результаты численного моделирования
Рассмотрим стационарное изотермическое течение разбавленного полимера в канале, имеющего препятствие в виде цилиндра, схема которого представлена на рис. 3.1. Рассмотрены два варианта течения: в первом рассматривается обтекание цилиндра в канале, стенки которого имеют скорость и (рис. 3.1, а), и второй, когда происходит пуазейлевское обтекание цилиндра в канале (рис. 3.1, б). На рисунке 2Н и 2К обозначают ширину канала и диаметр цилиндра, х,у — декартова система координат.
и
:и
2Ь
2ЬГ и ^^
а б
Рис. 3.1. Схемы расчетных областей при обтекании цилиндра в канале: а — однородный поток; б — течение Пуазейля
В центре расчетной области х = у = 0 расположен цилиндр, входная граница расположена на расстоянии от области Ь\ = 15Л,, а выходная на ¿2 = 30Н. В силу этого на входе задается полностью развитое течение полимера, определяющееся в безразмерном виде как
: и = 4 (4 - У2) , V = 0, Тхх = 2(1 - в) (|у) ,
Тху = (1 - вн
На выходе ставятся условия условия Неймана, предполагая, что здесь достигается установившееся течение
: ди = дгхх дх дх
дт
ху
дт
уу
0, р = 0.
дх дх
На твердых стенках канала ставятся условия прилипания жидкости и = = V = 0, а на оси ставится условие симметрии V = тху = 0. Моделирование течения проводилось при следующих значениях безразмерных параметров: Ев = 0, 01; в = 0,11; £ртт = 0,15; £ = 1, 0; ^е = 0, 01; 0, 5;1, 0;3, 0.
Реологическая конститутивная модель Фан-Тьен — Таннера
Одним из наиболее интересных и малоисследованных явлений реальных течений полимеров в канале при обтекании цилиндра является образование так называемого "отрицательного следа", возникающего только при движении груза в упруговязкой жидкости. Суть эффекта отрицательного следа заключается в том, что продольная скорость полимера за грузом восстанавливается от значения величины скорости груза до нулевого значения не сразу. В отличие от ньютоновского поведения, когда происходит монотонное снижение скорости, наблюдается снижение скорости до некоторой отрицательной величины с последующим ростом до нулевого значения. В случае обтекания цилиндрического груза за задней кромкой возникает эффект превышения скорости над скоростью обтекания цилиндра, что также характерно только для упруговязких жидкостей. В случае обтекания ньютоновской жидкостью указанный эффект отсутствует. Результаты наших расчетов значений продольной составляющей скорости для схем течения, представленных на рис. 3.1, приведены на рис. 3.2, из которого следует, что для схемы течения с подвижными стенками (рис. 3.2, а) наблюдается максимальное значение скорости при Шв = 0, 5, и оно составляет и = 1, 25, что означает 25 % превышение скорости над значением скорости на оси канала. При дальнейшем увеличении времени релаксации наблюдается некоторое уменьшение пика скорости, но процесс замедления до скорости на оси канала происходит на большем расстоянии.
для различных чисел Шв (1 — Шв = 0,01; 2 — 0,5; 3 — 1,0; 4 — 3,0): а — для схемы однородного течения (рис. 3.1, а); б — для течения Пуазейля (рис. 3.1, б)
Аналогичная ситуация наблюдается в поведении профиля скорости в следе за цилиндром при обтекании в канале потоком Пуазейля (рис. 3.2, б). Из этого же рисунка видно, что в случае, близком к ньютоновскому поведению при малых значениях времени релаксации ( Шв = 0, 01), ситуация
аналогична для обеих схем течения. Скорость монотонно восстанавливается из нулевого значения на цилиндре до скорости в основном потоке.
В дальнейшем определена такая важная характеристика течения упру-говязкой жидкости, как разность главных напряжений а\ — 02 = = л/+ 4тху, где N1 = тхх — туу. В соответствии с оптическим законом 01 — 02 = С(П1 —п2), где щ —П2 определяет разность коэффициентов преломления [4]. Таким образом, разность главных напряжений определяет оптическую неоднородность потока, что является следствием различной степени ориентации макромолекул полимерной жидкости. В свою очередь различная степень ориентации макромолекул в потоке существенным образом влияет на картину течения: а именно вытягивание молекулярных цепочек и их частичная ориентация в области течения вблизи твердой стенки способны привести к образованию дальнего порядка в расположении макромолекул. Полимерная жидкость, в этом случае, в пристенных слоях будет приобретать свойства упругого гуковского тела. Образующиеся за счет частичной ориентации участков макромолекул квазисшивки переводят расплав в пристенном слое из вязкотекучего в высокоэластическое состояние. При этом образующиеся вблизи стенки надмолекулярные структуры способны дополнительно уменьшить зазор между падающим грузом и стенкой контейнера, что неизбежно приведет к перераспределению структуры потока вблизи падающего груза.
На рис. 3.3 показано распределение величины разности главных напряжений в следе за цилиндром для обеих схем обтекания цилиндра, указанных на рис. 3.1. Из рис. 3.3 также следует, что при малых значениях числа Вайссенберга величина 01 —02 резко возрастает, а потом на расстоянии, равном 2 высотам канала от задней кромки цилиндра, уменьшается до нуля. С ростом значений числа Вайссенберга пик величины 01 — 02 снижается, однако релаксация главных напряжений происходит на гораздо большем расстоянии от задней кромки цилиндра (свыше 4 высот канала).
Реологическая конститутивная модель Олдройд-Б
Реологическая конститутивная модель упруговязкой жидкости характеризуется в первую очередь тем, что она предсказывает возникновение нормальных напряжений в сдвиговых течениях при постоянной сдвиговой вязкости. Таким образом, изначально постулируется, что полимерная жидкость, подчиняющаяся модели Олдройд-Б, не обладает свойством аномалией вязкости. На рис. 3.4 показано распределение продольной компоненты скорости в следе цилиндра. Для схемы течения Пуазейля в канале (рис. 3.1, б) численные расчеты, в соответствие с рис. 3.4, б, показывают монотонный характер восстановления скорости. В случае течения с подвижными стенками (рис. 3.1, а) данный эффект тоже имеет место, только выражен он очень незначительно при Шв = 0, 5. При увеличении числа Шв его проявление также исчезает. Таким образом, для обеих схем течения, представленных на рис. 3.1, эффект "отрицательного следа" незначителен
цилиндром для различных чисел Ше (1 — Ше = 0,01; 2 — 0,5; 3 — 1,0;
4 — 3,0): а — для схемы однородного потока (рис. 3.1, а); б — для течения
Пуазейля (рис. 3.1, б)
или практически отсутствует. Данная ситуация свидетельствует о том, что на возникновение эффекта "отрицательного следа" большее влияние оказывает наличие свойства аномалии вязкости упруговязкой жидкости.
Далее проанализируем распределение разности главных напряжений. Как было показано, в случае использования реологического соотношения РТТ, характеризующегося учетом аномалии вязкости, с ростом времени релаксации напряжений (ростом числа Ше), происходит некоторое уменьшение пика разности главных напряжений (рис. 3.3) с одновременным ростом расстояния, на котором а\ — ач релаксирует.
В случае использования реологического конститутивного соотношения Олдройд-Б ситуация меняется на противоположную. Как видно из рис. 3.5, с ростом времени релаксации напряжений в потоке (с ростом величины числа Вайссенберга) величина пика разности главных напряжений а\ — ач возрастает, а не убывает, как предсказывается моделью РТТ. Расстояние, на котором разность главных напряжений стремится к нулю, также возрастает. Данная ситуация свидетельствует о существенной роли свойства аномалии вязкости на распределение напряжений, а следовательно, и на общую картину течения вблизи цилиндра.
Рассмотрим картины течения, полученные для различных значений числа Вайссенберга. На рис. 3.6 обтекание цилиндра происходит слева направо. На рис. 3.6, а показаны линии уровня горизонтальной компоненты скорости и, на рис. 3.6, б показаны изобары, на рис. 3.6, в показаны линии уровня разности главных напряжений. Отметим, что, пользуясь возникновением эффекта двойного лучепреломления, вызванного течением полимерной жидкости вблизи цилиндра, можно экспериментально наблюдать изолинии п\ — пч, которые согласно оптическому закону пропорциональны раз-
а — для схемы однородного течения (рис. 3,1, а;) б — для течения Пуазейля (рис. 3.1, б)
Олдройд-Б для различных чисел Шв (1 — Шв = 0,01; 2 — 0,5; 3 — 1,0): а — для схемы однородного течения (рис. 3.1, а); б — для течения Пуазейля (рис. 3.1, б)
ности главных напряжений 01 — 02 [4]. Зона максимальной величины разности главных напряжений (рис. 3.6, в) расположена вдоль поверхности цилиндра, причем наибольшие их значения расположены на боковых поверхностях, где наблюдаются наибольшие градиенты скорости. На кривых распределения скорости потока и отчетливо проявляется зона превышения скорости (рис. 3.6, а). Также видно увеличение концентрации напряжения на осевой линии, а вследствие этого и 01 — 02 в зоне за цилиндром. Отметим, что максимальные значения безразмерных напряжений при этом уменьшаются, и их распределение по потоку становится более равномер-
ным. Как указывалось выше, реологическая модель упруговязкой жидкости Олдройд-Б применяется для полимерных жидкостей (растворов или расплавов), обладающих постоянной сдвиговой вязкостью. Поэтому в случае использования данной модели проявление эффекта "отрицательного следа" за цилиндром практически не наблюдается. Небольшое увеличение скорости было обнаружено только для схемы однородного течения (рис. 3.1, а) при Шв — 0, 5, однако при увеличении числа Шв этот эффект полностью исчезает.
Суммарная сила вязкого трения в проекции на продольную ось определяется по формуле
г2ж
■Рцгяс — /
ди \ (ди ду \
-р + 2^2дх + тхх) С08° + П2 ( дх + + Тху
К(Ш.
(3.1)
Однако большинство исследователей вычисляют коэффициент коррекции, определяемый по формуле
р .
С< тр _ тр ("< _ угяс
Сй • Г 8Ьссв — гvis, Са — , тт-
чпци
Коэффициент коррекции является достаточно важным интегральным параметром, показывающим влияние стенок канала и реологических свойств исследуемой жидкости на сопротивление жидкости падающему грузу. В табл. 3.1 приведены расчетные данные зависимости коэффициента коррекции от величины числа Вайссенберга.
Таблица 3.1
Значения коэффициента Са
о
'е Са Олдройд-Б Са РТТ Са Олдройд-Б Са РТ
рис. 3.1, а рис. 3.1, а рис. 3.1, б рис. 3.1, б
0,0 102,12 102,76 132,36 132,30
0,5 89,46 79,25 126,41 120,96
1,0 81,96 70,45 121,12 118,77
3,0 — 61,26 — 105,65
Из табл. 3.1 следует, что для заданных геометрических соотношений интегральное влияние неньютоновских свойств приводит к уменьшению величины коэффициента коррекции. В результате расчетов коэффициента коррекции при обтекании цилиндров установлено, что свойство аномалии вязкости играет более значительную роль в снижении сопротивления цилиндра. Уменьшение сопротивления цилиндра по модели РТТ, учитывающей данное свойство, является более значительным. Коэффициент сопротивления для плоского обтекания, когда поток ограничен двумя неподвижными стенками (схема на рис. 3.1, б), превышает величину сопротивления в случае относительного движения стенок канала (схема течения на рис. 3.1, а) приблизительно на 20 %.
Рис. 3.6. Картина течения жидкости РТТ для рис. 3.1, а при Ш в = 1 в виде:
а — изолиний продольной скорости; б — изобар; в — изолиний о\ — о^
Заключение
Результаты численного моделирования стационарного обтекания цилиндра в канале с неподвижными стенками и в канале, стенки которого движутся с некоторой скоростью, показали, что в обоих случаях вязкоупругая жидкость в следе за цилиндром проявляет нелинейное поведение, отличное от стоксового обтекания ньютоновской жидкостью. Данное свойство проявляется в образовании так называемого "отрицательного следа", профиль скорости за цилиндром носит немонотонный характер. Проведенное численное исследование данного эффекта показывает рост зоны "отрицательного следа" при увеличении времени релаксации жидкости, что позволяет сделать предположение, что данный эффект связан с проявлением эффекта нормальных напряжений, а также аномалией вязкости вязкоупругого материала в сдвиговом течении в такой форме. Расчеты для модели жидкости Олдройд-Б, обнаруживают очень слабое проявление данного эффекта, что подтверждает утверждение о значительном влиянии в этом обтекании свойства аномалии вязкости вязкоупругого материала.
Отметим также, что в данной работе предлагается новый численный алгоритм, использованный для моделирования вязкоупругих течений полимеров в каналах различной формы. Представлены результаты расчетов вязко-упругого течения полимера в канале при обтекании цилиндра. Рассмотрен такой неньютоновский эффект, как образование отрицательного следа при обтекании цилиндра в канале. Численные результаты хорошо согласуются с физическими представлениями о течениях аномально вязких жидкостей в сложных каналах. Полученные результаты позволяют утверждать, что, так как в данном эффекте влияние учета аномалии вязкости является очень значительным, если не определяющим, выбор реологической модели, позволяющей учитывать аномалию вязкости для изучения данного эффекта, является необходимым фактором.
Литература
[1] Хаппель Дж., Бренер Г. Гидродинамика при малых числах Рейнольд-са. М.: Мир, 1976. 623 с.
[2] Dou H.S., Phan-Thien N. Negative wake in the uniform flow past a cylinder // Rheol. Acta. 2003. V. 42. P. 383-409.
[3] Phan-Thien N., Tanner R.I. A new constitutive equation derived from network theory // J. Non-Newtonian Fluid. Mech. 1979. V. 2. P. 353-365.
[4] Bird R.B., Armstrong R.C., Hassager O. Dynamics of Polymeric Liquids. vol.1. Fluid Mechanics. 2nd ed. New York: John Wiley and Sons, 1987. 565 с.
[5] Назмеев Ю.Г. Гидродинамика и теплообмен закрученных потоков реологически сложных жидкостей. М.: Энергоатомиздат, 1996. 304 c.
[6] Коннор Дж., Бреббиа К. Метод конечных элементов в механике жидкости. Л.: Судостроение, 1979. 264 с.
[7] High-order finite volume schemes for viscoelastic flows / B. A. Snigerev [et al.] // J. Comput. Phys. 2004. V. 199. P. 16-40.
[8] Luo X.L. Operator splitting algorithm for viscoelastic flow and numerical analysis for the flow around sphere in tube // J.Non-Newtonian Fluid. Mech. 1996. V. 48. P. 57-75.
[6] Писсанецки В.И. Технология разреженных матриц. М.: Мир, 1988. 412 с.
[9] Arigo M.T., McKinley G.H. An experimental investigation of negative wake behind spheres settling in a shear-thinning visoelsctic fluid // Rheol. Acta. 1998. V. 37. P. 307-327.
Поступила в редакцию 24/VT/2009; в окончательном варианте — 24/VT/2009.
THE INVESTIGATION OF NON-NEWTONIAN EFFECT AT FLOW PAST CYLINDER BY THE CRALLING FLOW OF VISCOELASTIC FLUID IN THE FLAT CHANNEL WITH FALLING WEIGHT
© 2009 B.A. Snigerev2
The flow of polymeric liquid in the flat channel with falling weight is considered. The rheological constitutive models for fluids are the following: the viscoelastic model of Oldroyd-B fluid, relaxation model of Phan-Thien-Tanner, which prognosticate polymer solution viscosity anomaly and presence of parallel viscosity in flows, conditioned by normal stress. For this models the numerical results of the task about the falling weight in viscometer are received at different boundary conditions on the walls of the channel. The results have good agreements with numerical and experimental data of other authors.
Key words: viscoelastic flow, flat channel, finite element method.
Paper received 24/VI/2009. Paper accepted 24/VI/2009.
2Snigerev Boris Alexandrovich (snigerev@mail.knc.ru), Lab. of Modelling of Technological Processes, Institute of Mechanics and Engineering, RAS, Kazan, 420111, Russia.
