CC BY
60
раздел МАТЕМАТИКА и МЕХАНИКА
УДК 532.593
УДАР СТРУИ ПО ТОНКОМУ СЛОЮ ЖИДКОСТИ НА СТЕНКЕ © А. А. Аганин, Т. С. Гусева*
Институт механики и машиностроения Казанского научного центра РАН Россия, Республика Татарстан, 420111 г. Казань, ул. Лобачевского, 2/31.
Тел.: +7 (843) 231 91 14.
*Етай: [email protected]
Изучается удар высокоскоростной струи воды по тонкому водяному слою (прослойке) на плоской жесткой стенке. Параметры задачи характерны для ударного воздействия на стенку струи, возникающей на поверхности кавитационного пузырька при его схлопывании вблизи твердой стенки. Для прослойки толщиной 0.06 радиуса струи проводится анализ влияния скорости струи в диапазоне от 150 до 350 м/c. При таких скоростях удара демпфирующее влияние прослойки относительно мало, и на стенке возникают давления того же порядка, что и в случае без прослойки. Исследование проводится с применением прямого численного моделирования с описанием движения жидкости и газа уравнениями газовой динамики. Показано, что при скорости струи больше 200 м/с распределение давления на стенке существенно неоднородно, как и в случае без прослойки. При этом максимум давления также реализуется на краю области импульсного воздействия и превышает давление гидроудара.
Ключевые слова: удар струи, прослойка жидкости на стенке, импульс давления на стенке.
Введение
Задача об ударе жидкой массы (струи или капли) по тонкому слою жидкости на стенке актуальна для многих приложений. Такие удары реализуются при работе паровых турбин, при полете летательных аппаратов в условиях дождя, при коллапсе кавитационных пузырьков вблизи стенки [1]. Удар струй, образующихся при схлопывании кавитацион-ных пузырьков, считается одной из возможных причин широко известного явления разрушения поверхностей тел, работающих в условиях кавитации [1-3].
Экспериментальное исследование удара по стенке струи, образующейся при схлопывании пузырька (рис. 1), осложнено малостью пространственных и временных масштабов явления. Так, радиус конца струи при ударе Я ~ 1-10 мкм, время прохождения волны разрежения от периферии струи до ее оси х = Я / С ~ 1-10 нс [1], С - скорость звука в жидкости. Еще более сложной задачей для экспериментальных исследований является удар струи по тонкой прослойке жидкости на стенке, когда пузырек схлопывается на небольшом удалении от стенки. При малой толщине прослойки ее демпфирующее влияние достаточно мало, и уровни давления на стенке могут быть сравнимыми с тем, что реализуются при ударе по сухой стенке.
Характерными особенностями импульсного воздействия струи на стенку являются большие и
быстрые деформации границы раздела между жидкостью и газом, изменение ее связности, ударные волны в жидкости и газе. Некоторые из этих особенностей иллюстрирует рис. 2, где приведены фотографии экспериментального моделирования удара струи по стенке [5]. В результате удара со скоростью 150 м/с в моделирующем струю геле сначала возникает ударная волна (рис. 2Ь), а затем реализуется боковое растекание геля по поверхности образца (рис. 2с).
Рис. 1. Фотография схлопывающегося у стенки пузырька с бьющей по стенке струей [4].
Ударные волны и большие деформации контактной границы жидкость-газ накладывают довольно жесткие ограничения на выбор математической модели и метода расчета удара струи по стенке. Так, при выборе модели нужно учитывать сжимаемость жидкости. При выборе численного метода нужно учитывать, что применение лагранжева подхода будет осложнено большими деформациями частиц среды, особенно на стадии растекания струи.
Рис. 2. Удар летящего вверх металлического образца по неподвижному слою геля дискообразной формы [5].
Применение смешанного лагранжево -эйлерова подхода с возможностью перестройки сетки как, например, в [6], позволяет несколько ослабить эту проблему. Однако по мере увеличения скорости и величины деформации межфазной границы эффективность перестройки снижается. Использование неструктурированных перестраивающихся сеток, как, например, в [7], значительно усложняет алгоритм.
В настоящей работе изучается удар высокоскоростной струи жидкости по тонкой жидкой прослойке, находящейся на жесткой стенке. Параметры задачи характерны для ударного воздействия на стенку струи, возникающей на поверхности кавита-ционного пузырька при его схлопывании вблизи твердой стенки. Рассматривается диапазон скоростей удара, при которых демпфирующее влияние прослойки относительно мало, и на стенке возникают давления того же порядка, что и в случае без прослойки. Используется методика расчета, описанная в [8-11], основанная на прямом численном моделировании с описанием движения жидкости и газа уравнениями газовой динамики. Применяется метод CIP-CUP (Constrained Interpolation Profile - Combined Unified Procedure) [12] в сочетании с подходом без явного выделения межфазной границы [13]. Метод CIP-CUP достаточно эффективен при умеренной интенсивности ударных волн и контактных границах с сильным перепадом акустического импеданса (в том числе и границы типа газ-жидкость) [13]. Применяется неструктурированная динамически-адаптивная soroban-сетка [14].
1. Постановка задачи
Рассматривается осесимметричная задача об ударе цилиндрической струи жидкости с полусферическим концом по тонкому равномерному слою жидкости (прослойке) на плоской жесткой стенке (рис. 3a). Радиус струи R, толщина прослойки d, жидкость в струе и прослойке одна и та же (вода). Струя ударяет по прослойке ортогонально поверхности стенки.
150 200 250 300 350 Рис. 3. Струя и прослойка жидкости на стенке в начале удара (а) и изменение давления гидроудара pwh в рассматриваемом диапазоне скорости струи V (b).
Исследуется влияние скорости струи V в диапазоне 150 < V < 350 м/с при толщине прослойки d / R = 0.06. Такие скорости струи и толщина прослойки характерны для ударного воздействия на стенку кавитационного пузырька. Близкие скорости струй получены при схлопывании кавитационного
пузырька в [16]. Они соответствуют диапазонам скоростей и толщин, при которых на стенке реализуются давления порядка давления гидроудара (water hammer pressure) pwh = ploDV, где plo - невозмущенная плотность жидкости, D - скорость ударной волны, распространяющейся вверх по струе (для воды D и С + 2V) [1, 15]. При непосредственном ударе струи по сухой стенке (без прослойки жидкости) давление pwh возникает на стенке в начальный момент удара, а затем в большей части области воздействия реализуются давления, близкие к pwh. Зависимость давления гидроудара от скорости струи для водяной струи при pL0 =103 кг/м3, С =1500 м/с приведена на рис. 3b. В рассматриваемом диапазоне скоростей давление гидроудара достигает довольно больших значений (максимум составляет около одного гигапаскаля), что представляет значительный интерес для приложений.
2. Основные положения метода расчета
В методике расчета применяется подход без явного выделения межфазной границы. Вводится функция-идентификатор ф, такая, что в области жидкости ф =1, в области газа ф = 0, а в малой окрестности их контактной границы ф непрерывна и монотонно меняется от 0 до 1. Контактная граница при этом рассматривается как переходная область, занятая смешанной средой. Изменение идентификатора описывается уравнением переноса
ф1 + u-Уф = 0 (1)
Для сквозного расчета динамики жидкости и газа используется метод CIP-CUP [12]. Он базируется на уравнениях движения сжимаемой среды в неконсервативных переменных (плотность p, скорость u, давление p), что и позволяет осуществлять сквозной расчет без явного выделения межфазной границы [17]. Динамика контактирующих жидкости и газа без учета эффектов вязкости и теплопроводности описывается уравнениями, которые в векторном виде можно записать как
ft + u Vf = G. (2)
Здесь f = (p, u, p)T, G = (-pV-u, -p-1Vp, -pCs2V-u)T, где p - плотность, u - скорость, p - давление, Cs = фСя1+(1 - ф) Cs2 - скорость звука, CSi = [r,(p+BI)/p]1/2, i = 1 соответствует жидкости, i = 2 - газу, Ti, B1 - константы уравнения состояния жидкости в форме Тэта, Г2 = у, B2 = 0, у - показатель адиабаты газа.
Численное решение системы (2) разбивается на конвективную и неконвективную стадии
f * - f"
- + u-Vf = 0,
At" f"+l _f*
At"
= G.
(3)
(4)
Для расчета группы уравнений переноса (3) применяется полулагранжев метод CIP (Constrained
Interpolation Profile) [13]. Решение уравнения переноса для функции f в узле x в момент tn+1 = tn + Atn можно приближенно записать в виде f (x-uAtn), где аргумент представляет собой отправную точку - положение «лагранжевой частицы», которая за время Atn перемещается в узел x. Значение переносимой величины f в отправной точке определяется CIP-интерполяцией с использованием известных узловых значений f и ее пространственных производных fx, fy, fxy в момент tn.
В рамках метода UP (Unified Procedure) [12], неконвективная часть (4) сводится к следующему уравнению на давление
и+1 *
* * ,n2
Р cs At
= V-
rVpn+1 ^
V-u
' Atn
(5)
которое решается итерационным методом. Затем пересчитываются значения скорости un+1 и плотности pn+1. При решении задач с интенсивными ударными волнами метод CIP-CUP применяется в сочетании с искусственной вязкостью.
Рис. 4. Двумерная soroban-сетка, j - номер направляющей линии, i - номер узла.
Расчеты проводятся с использованием адаптивной soroban-сетки (от японского soroban - счеты) [14]. В двумерном случае soroban-сетка - это набор узлов, расположенных на параллельных направляющих линиях (рис. 4). На каждом временном шаге строится новая сетка, адаптирующаяся к решению и никак не связанная со старой сеткой. При этом взаимное расположение и количество как направляющих, так и узлов на них может изменяться.
При построении сетки сначала рассчитывается монитор-функция
1/2
M(x, y,t) =(l + а( fx2 + fy2 )) + p(| fxJ\ + fyy\ )•
Здесь f - какой-либо параметр решения, a, p -положительные константы. Направляющие и узлы сетки сгущаются в зонах с большими значениями монитор-функции и разрежаются там, где ее значения малы. При построении конечно-разностных аппроксимаций производных на soroban-сетке требуемые значения параметров в узлах шаблона схемы определяются CIP-интерполяцией [14].
Разработанные на основе этой методики алгоритмы и программы были протестированы на многих задачах в одномерной и двумерной постановке. Численные решения сравнивались с известными
аналитическими решениями и численными решениями других авторов [10, 11].
3. Результаты
Для изучения влияния скорости струи рассчитывались 5 вариантов удара: V = 150, 200, 250, 300 и 350 м/с. Основные данные получены на достаточно мелкой сетке М0 с минимальным размером ячеек 0.004 ё. Для уточнения отдельных деталей рассматриваемого процесса применялись сетки М1 и М2 с минимальными размерами ячеек, меньшими, чем у сетки М0, в 2 и 4 раза, соответственно.
На рис. 5 подробно представлен случай V = 350 м/с. При ударе струи по слою жидкости возникают две ударные волны: первая уходит вверх по струе, вторая - вниз по прослойке. По достижении второй ударной волной стенки возникает отраженная от стенки ударная волна (рис. 5а). На начальном этапе края первой и второй ударных волн примыкают к кромке контакта струи и прослойки (рис. 5а, точка А). Это обусловлено высокой скоростью расхождения кромки, поскольку в ее окрестности угол наклона поверхности струи к поверхности прослойки еще мал. Видно, что этот угол будет постепенно увеличиваться, поэтому скорость расхождения кромки будет снижаться. В момент 3 (рис. 5Ь) ударные волны уже обгоняют кромку. Примерно до момента 3 реализуется регулярное отражение падающей на стенку ударной волны.
К моменту 4 (рис. 5с) взаимодействие ударных волн с межфазной границей вблизи кромки области контакта струи и прослойки приводит к развитию очень тонкого кольцевого выплеска жидкости. Примерно в это же время отражение падающей на стенку ударной волны, которая начинает перемещаться в радиальном направлении по прослойке, переходит в нерегулярный режим (с образованием ножки Маха).
К моменту 5 (рис. 5ё) кольцевой выплеск жидкости между струей и прослойкой, а также ножка Маха в окрестности стенки становятся выраженными. Отраженная от стенки ударная волна достигает основания кольцевого выплеска жидкости на кромке области контакта струи и прослойки. В результате этого одна часть фронта отраженной ударной волны уходит в прослойку, а другая часть - в струю.
К моменту 6 (рис. 5е) фронт отраженной ударной волны догоняет фронт первой ударной волны вблизи боковой поверхности струи. В области прослойки реализуется взаимодействие ударной волны со свободной поверхностью с образованием волны разрежения, распространяющейся к стенке. В результате ее отражения от стенки в прослойке возникает область с довольно большими отрицательными давлениями (момент 7, рис. 5 f ), что на практике означало бы возникновение кавитационных полостей. Используемая в настоящей работе математическая модель не позволяет описать динамику таких полостей и оценить их последствия. Но, опираясь на
"1—'—I—1—I—»—I—1—I
О 0.1 0.2 0.3 0.4 //т 0.5
Рис. 5. Удар струи со скоростью V = 350 м/с: (a-f) - поля давления (красная кривая - межфазная граница) в половине осевого сечения; (g) - радиальные профили давления на стенке (кривые 2-7 соответствуют (a-f )); (h) - временные зависимости максимального давления на стенке (кривые M0-M2) и в центре области воздействия (кривая pc), моменты 1-7
отмечены точками, т = R/C = 6.8 нс, pwh = 0.77 ГПа.
обширные исследования кавитационного разрушения, можно предположить, что это повышает разрушительный потенциал ударного воздействия струи.
Импульсное воздействие на стенку начинается, когда ее достигает ударная волна, расходящаяся от кромки контакта струи и прослойки. В начале воздействия давление в малой окрестности центра области отражения (области ударного воздействия) увеличивается скачком до 0.75 ркк (рис. 5g, к). В последующем давление на стенке в центре области воздействия монотонно уменьшается, а на периферии оно сначала возрастает (примерно до момента 3), а затем убывает. Известно, что аналогичная неоднородность наблюдается в профилях давления на стенке при отсутствии прослойки. После момента 4 на периферии области воздействия начинает развиваться обусловленный волной разрежения локальный минимум давления, который может достигать достаточно больших отрицательных величин. Это возможно лишь при наличии прослойки.
Представленные на рис. 5а^ результаты получены на сетке М0. Применение более мелких сеток
М1 и М2 позволило несколько уточнить численное решение. В профилях давления на стенке это уточнение проявляется лишь как увеличение высоты периферийных пиков давления в небольшой кольцевой окрестности r/R и 0.25 в коротком интервале 0.12< t/т < 0.18. В остальном картина импульсного нагружения в основной части области воздействия, в том числе в ее центре, при уменьшении минимального шага сетки не изменяется. Как видно (рис. 5h), даже столь мелкие сетки, как М1 и М2, не позволяют окончательно оценить величину пиковых давлений на стенке в силу их крайне малых пространственно -временных масштабов.
Уменьшение скорости струи в рассматриваемом диапазоне 150 < V < 350 м/с не приводит к существенному изменению качественных особенностей динамики жидкости в струе и прослойке. Возникает ряд количественных изменений. С уменьшением V снижаются уровни давления на кромке контакта струи и прослойки, интенсивность расходящихся от кромки ударных волн. В результате пони-
О 0.2 0.4 r/R 0.6 0 0.1 0.2 0.3 tlx 0.4
Рис. 6. То же, что и на рис. 5 g, h, но для удара струи со скоростью V=150 м/с, pwh = 0.27 ГПа.
жаются уровни давления на стенке, и изменяется характер его распределения. В частности, при V = 150 м/с (рис. 6) максимум давления на стенке примерно равен 0.6 р„н. При этом распределение давления на стенке в начале воздействия (до формирования локального минимума в окрестности периферии, моменты 4, 5) остается близким к однородному. При этом, как видно на рис. 6Ь, кривые максимального давления на стенке, полученные на сетках М0-М2, практически совпадают.
— P'Pvh
-
..»--------■•
- ............
- -М2 .....МО
-------Ml 1 1 .........Рс 1 1 V, м/с 1 1 1
150 200 250 300 350
Рис. 7. Зависимости максимумов давления на стенке во всей области воздействия струи (кривые M0-M2) и в ее центре (кривая pc) от скорости струи V.
Из рис. 7 следует, что при ударе водяной струи по водяной прослойке толщиной d = 0.06R максимальные давления на стенке при V > 300 м/c достигают уровня давления гидроудара pwh (характерный уровень давлений при ударе по сухой стенке). Расхождение кривых максимального давления во всей области воздействия и в ее центре с ростом V обусловлено возникновением на периферии этой области острых пиков давления и увеличением их высоты.
Заключение
Проведено исследование удара высокоскоростной водяной струи по тонкому водяному слою (прослойке) на жесткой стенке. Толщина прослойки составляет 0.06 от радиуса струи, скорость струи V варьируется в диапазоне от 150 до 350 м/с. Подобные параметры струи и прослойки типичны для ударного воздействия на стенку струи, возникающей на поверхности кавитационного пузырька при его схлопывании вблизи твердой стенки. Установлено следующее.
1. При V > 300 м/c максимальные давления на стенке достигают уровня давления гидроудара Pwh = PlqDV (pL0 - невозмущенная плотность жидкости, D - скорость ударной волны, распространяющейся вверх по струе).
2. Максимум давления в центре области ударного воздействия на стенку меняется от 0.6 pwh при 150 м/c до 0.75 pwh при 350 м/c. Он достигается в начале отражения ударной волны от стенки. В последующем давление в центре области ударного воздействия плавно уменьшается.
3. При V < 200 м/c распределение давления на стенке в области ударного воздействия (до начала взаимодействия ударной волны с поверхностью прослойки) близко к однородному, а с ростом V становится все более неоднородным. На периферии области воздействия возникают острые пики давления, высота которых по мере роста V увеличивается.
4. В результате взаимодействия распространяющейся по прослойке ударной волны с поверхностью прослойки образуется волна разрежения. При ее отражении от стенки в прослойке возникает область с довольно большими отрицательными давлениями. На практике это может приводить к возникновению кавитации, которая, как известно, повышает разрушительный потенциал воздействия жидкости на тело.
Работа выполнена при поддержке Российского фонда фундаментальных исследований (проект № 16-01-00433).
ЛИТЕРАТУРА
1. Bourne N. K. On impacting liquid jets and drops onto polymethylmethacrylate targets // Proc. R. Soc. A. 2005. Vol. 461. Pp. 1129-1145.
2. Kornfeld M., Suvorov L. On the destructive action of cavitation // J. Appl. Phys. 1944. Vol. 15. Pp. 495-506.
3. Hawker N. A., Ventikos Y. Interaction of a strong shockwave with a gas bubble in a liquid medium: a numerical study // J. Fluid Mech. 2012. Vol. 701. Pp. 59-97.
4. Crum L. A. Surface oscillations and jet development in pulsating bubbles // J. Physique. 1979. Vol. 40. No. 11. Pp. C8-285-288.
5. Dear J. P., Field J. E. High-speed photography of surface geometry effects in liquid/solid impact // J. Appl. Phys. 1988. Vol. 63. Pp. 1015-1021.
6. Аганин А. А., Халитова Т. Ф., Хисматуллина Н. А. Метод численного решения задач сильного сжатия несферического кавитационного пузырька // Вычислительные технологии. 2010. Т. 15. №1. C. 14-32.
7. Головачeв Ю. П., Ноткина Е. А., Чижов А. В., Шмидт А. А. Расчет ударно-волновых течений со свободными поверхностями // ЖВММФ. 2001. Т. 41. №1. C. 157-167.
8. Аганин А. А., Гусева Т. С. Численное моделирование контактного взаимодействия сжимаемых сред на эйлеровых сетках // Ученые записки Казанского университета. Сер. Физ.-матем. науки. 2012. Т. 154. Кн. 4. С.74-99.
9. Аганин А. А., Гусева Т. С. Расчет контактного взаимодей- 13. ствия сжимаемых сред без явного выделения межфазных границ // Вестник Башкирского Университета. 2013. Т. 18. №3. С. 646-661. 14.
10. Аганин А. А., Гусева Т. С. Численное моделирование динамики неоднородных сжимаемых сред на основе метода CIP-CUP на адаптивных Soroban-сетках // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки 2014. Т. 156. Кн. 2. С. 55-72. 15.
11. Аганин А. А., Гусева Т. С. Методика расчета волн в жидкости и газе методом CIP-CUP с применением динамически-адаптивных Soroban-сеток // Вестник Башкирского Уни- 16. верситета. 2014. Т.19. №2. С. 368-380.
12. Yabe T., Wang P. Y. Unified numerical procedure for compressible and incompressible fluid // J. Phys. Soc. Japan. 1991. 17. Vol. 60. No. 7. Pp. 2105-2108.
Yabe T., Xiao F., Utsumi T. The constrained interpolation profile method for multiphase analysis // J. Comput. Phys. 2001. Vol. 169. No. 2. Pp. 556-593.
Takizawa K., Yabe T., Tsugawa Y., Tezduyar T. E., Mizoe H. Computation of free-surface flows and fluid-object interactions with the CIP method based on adaptive meshless Soroban grids // Comput. Mech. 2007. Vol. 40. Pp. 167-183. Heymann F. J. High-speed impact between a liquid drop and a solid surface // J. Appl. Phys. 1969. Vol. 40. No. 13. Pp. 51135122.
Воинов О. В., Воинов В. В. О схеме захлопывания кавита-ционного пузырька у стенки и образовании кумулятивной струйки // Докл. АН СССР. 1976. Т. 227. №1. С. 63-66. Abgrall R., Karni S. Computations of compressible multifluids // J.Comp.Phys. 2001. Vol. 169. Pp. 594-623.
Поступила в редакцию 25.04.2016 г.
COLLISION OF HIGH-SPEED WATER JET WITH A THIN LIQUID LAYER ON A RIGID WALL
© A. A. Aganin, T. S. Guseva*
Institute of Mechanics and Engineering, Kazan Science Center of RAS 2/31 Lobachevsky St., 420111 Kazan, Republic of Tatarstan, Russia.
Phone: +7 (843) 231 91 14.
*Email: [email protected]
A mechanics of collision of a high-speed water jet with a thin water layer on a plane rigid wall was studied. The problem parameters are typical for the water jet arising on the surface of a cavitation bubble during its collapse near a rigid wall. The analysis of the collision of the water jet having velocity V with the layer of thickness equal to 0.06 times of the jet radius was carried out. The range of the impact velocity from 150 to 350 m/s was considered. The was conducted by applying the direct numerical simulation with describing the motion of the liquid and the gas by the gas dynamics equations. It has been shown that the pressure maximum in the center of the area of collision on the wall grows from 0.6 Pwh at V = 150 m/s to 0.75 Pwh at V = 350 m/s (Pwh is the water-hammer pressure). For V< 200 m/s the pressure distribution on the wall (until the beginning of interaction of the shock wave with the surface of the liquid layer) is close to uniform, while it becomes more and more nonuniform with the growth in V. Sharp spikes appear on the periphery of the impact area, the height of which increases with increasing V. For the jet velocity greater than 300 m/s the maximum pressures on the wall are higher than the Pwh. As a result of interaction of the shock wave propagating in the liquid layer with its surface, a rarefaction wave is formed. During its reflection from the wall a region with large negative pressures arises. In practice it can lead to appearance of cavitation, which is known to increase the damaging effect of the water on the solid body.
Keywords: jet impact, liquid layer, pressure pulse on wall, shock wave in jet.
Published in Russian. Do not hesitate to contact us at [email protected] if you need translation of the article.
REFERENCES
1. Bourne N. K. Proc. R. Soc. A. 2005. Vol. 461. Pp. 1129-1145.
2. Kornfeld M., Suvorov L. J. Appl. Phys. 1944. Vol. 15. Pp. 495-506.
3. Hawker N. A., Ventikos Y. J. Fluid Mech. 2012. Vol. 701. Pp. 59-97.
4. Crum L. A. J. Physique. 1979. Vol. 40. No. 11. Pp. C8-285-288.
5. Dear J. P., Field J. E. J. Appl. Phys. 1988. Vol. 63. Pp. 1015-1021.
6. Aganin A. A., Khalitova T. F., Khismatullina N. A. Vychislitel'nye tekhnologii. 2010. Vol. 15. No. 1. Pp. 14-32.
7. Golovachev Yu. P., Notkina E. A., Chizhov A. V., Shmidt A. A. ZhVMMF. 2001. Vol. 41. No. 1. Pp. 157-167.
8. Aganin A. A., Guseva T. S. Uchenye zapiski Kazanskogo universiteta. Ser. Fiz.-matem. nauki. 2012. Vol. 154. Kn. 4. Pp. 74-99.
9. Aganin A. A., Guseva T. S. Vestnik Bashkirskogo Universiteta. 2013. Vol. 18. No. 3. Pp. 646-661.
10. Aganin A. A., Guseva T. S. Uchen. zap. Kazan. un-ta. Ser. Fiz.-matem. nauki 2014. Vol. 156. Kn. 2. Pp. 55-72.
11. Aganin A. A., Guseva T. S. Vestnik Bashkirskogo Universiteta. 2014. Vol. 19. No. 2. Pp. 368-380.
12. Yabe T., Wang P. Y. J. Phys. Soc. Japan. 1991. Vol. 60. No. 7. Pp. 2105-2108.
13. Yabe T., Xiao F., Utsumi T. J. Comput. Phys. 2001. Vol. 169. No. 2. Pp. 556-593.
14. Takizawa K., Yabe T., Tsugawa Y., Tezduyar T. E., Mizoe H. Comput. Mech. 2007. Vol. 40. Pp. 167-183.
15. Heymann F. J. J. Appl. Phys. 1969. Vol. 40. No. 13. Pp. 5113-5122.
16. Voinov O. V., Voinov V. V. Dokl. AN SSSR. 1976. Vol. 227. No. 1. Pp. 63-66.
17. Abgrall R., Karni S. J.Comp.Phys. 2001. Vol. 169. Pp. 594-623.
Received 25.04.2016.
