Расчет ударных волн в переменных плотность, скорость, давление Текст научной статьи по специальности «Физика»

А. А. Аганин, Т. С. Гусева

РАСЧЕТ УДАРНЫХ ВОЛН В ПЕРЕМЕННЫХ ПЛОТНОСТЬ, СКОРОСТЬ, ДАВЛЕНИЕ

Ключевые слова: эйлерова сетка, метод CIP-CUP, ударные волны, искусственная вязкость.

Изучены особенности расчета ударных волн в переменных плотность, скорость, давление с помощью разработанной авторами методики на основе метода CIP-CUP. Одно из важных достоинств этого метода - возможность сквозного расчета ударных волн и контактных границ типа газ-жидкость. Показано, что при расчете ударных волн в жидкости с существенным скачком энтропии возникают нефизичные осцилляции в окрестности фронта, а характеристики потока за фронтом и скорость его распространения определяются с большими погрешностями. Введение искусственной вязкости позволяет устранить все эти проблемы. Если искусственная вязкость вводится в уравнения для скорости и давления, то в численном решении на фронте ударной волны помимо массы и импульса сохраняется полная энергия, а если лишь в уравнение для скорости -то энтропия.

Keywords: Eulerian grid, CIP-CUP method, shock waves, artificial viscosity.

Features of computing shock waves in variables density, velocity and pressure by a numerical technique developed by the authors based on the method CIP-CUP have been investigated. One of important merits of the method is that it allows one to evaluate, without explicit separation, shock waves and contact discontinuities of gas-liquid type. It has been shown that computation of shock waves in liquid with essential entropy jump leads to unphysical oscillations in the vicinity of their fronts and gives significant errors in characteristics of the flow behind the shock front and velocity of its propagation Introduction of artificial viscosity allows one to remove all these problems. Introducing artificial viscosity in the equations for velocity and pressure gives numerical solution with total energy conservation (apart from conserving mass and momentum) on the shock front, while introducing it in the velocity equation only leads to solution with entropy conservation.

Введение

Многие технологические процессы представляют собой течения сжимаемых сред с контактными границами, которые претерпевают существенные деформации, фрагментацию, объединение. В качестве примера можно привести взаимодействие капель с пульсирующим газовым потоком в цилиндрической трубе в процессе горения жидкого топлива [1], диспергирование пузырей и последующее взаимодействие сплошной и дисперсной фаз в пульсационных аппаратах проточного типа, служащих для интенсификации процессов тепло- и массопереноса [2]. В некоторых процессах большая деформация контактных границ сопровождается интенсивными ударными волнами. Это наблюдается, например, при схлопывании кавитационных пузырьков с формированием кумулятивных струй вблизи стенок, в частности, при ультразвуковой очистке поверхностей [3], при сильном сжатии кавитационных пузырьков в таких технологиях как получение применяемой при производстве шин полимерной серы [4], гидродинамическая кавитационная обработка воды с целью ее обеззараживания [5].

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

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

Если контактные границы сильно деформируются, распадаются на части, или, наоборот, возникают в результате объединения отдельных структур, то их явное выделение становится весьма затруднительным. В таких случаях контактные границы удобнее рассчитывать на неподвижных (эйлеровых) сетках сквозным образом (аналогично сквозному счету ударных волн). При этом для неявного отслеживания контактных границ применяется какая-либо пространственная функция-идентификатор среды, принимающая в областях, занятых разными средами, постоянные различающиеся значения, а в окрестности контактной границы непрерывно и монотонно изменяющаяся между этими значениями. Переходная область в окрестности контактной границы с промежуточными значениями функции-идентификатора может интерпретироваться как некая смесь контактирующих сред, физический смысл которой может быть и не определен. Изменение функции-идентификатора описывается уравнением переноса. Положение контактной границы можно

устанавливать, например, по линиям определенного уровня или области максимального градиента функции-идентификатора.

Для расчета задач с большими деформациями контактных границ и ударными волнами, в явном выделении которых нет необходимости, наиболее естественным является подход, в котором и ударные волны, и контактные границы рассчитываются на эйлеровой сетке сквозным образом (т.е. отслеживаются неявно). Среди методов сквозного счета ударных волн с их размыванием с помощью схемной вязкости одним из наиболее популярных в настоящее время является метод С.К. Годунова (и его многочисленные модификации). Он широко применяется для численного решения задач динамики газа, жидкости, упруго-пластических тел и т.д. В методе С.К. Годунова исходные уравнения обычно записываются в дивергентной форме относительно массы, импульса и полной энергии единицы объема. Их разностные аппроксимации выражают законы сохранения массы, импульса и полной энергии отдельных ячеек расчетной сетки. В задачах контактного взаимодействия при явном выделении контактных границ метод С. К. Годунова применяется естественным образом, без внесения в его алгоритм каких-либо существенных изменений. Однако сквозной расчет без явного выделения контактных границ с использованием этого метода оказывается проблематичным. Согласно [6] использование консервативных схем сквозного счета приводит к появлению нефизичных осцилляций давления большой амплитуды в окрестности контактных границ.

Избежать таких осцилляций можно, воспользовавшись неконсервативными схемами сквозного счета в переменных плотность, скорость, давление [6]. Выполнение соотношений Гюгонио на ударных волнах при использовании неконсервативных схем сквозного счета обеспечивается с помощью искусственной вязкости.

Сравнительно недавно был предложен метод CIP-CUP (Constrained Interpolation Profile - Combined Unified Procedure) [7], позволяющий сквозным образом на эйлеровых сетках рассчитывать ударные волны и контактные границы с сильным перепадом акустического импеданса. В работе [8] показано, что этот метод является эффективным для расчета задач с большими деформациями контактных границ и интенсивными ударными волнами. В настоящей работе исследуются особенности расчета ударных волн в переменных плотность, скорость, давление с помощью разработанной авторами методики расчета на основе метода CIP-CUP.

1. Уравнения динамики сжимаемой среды

Уравнения динамики сжимаемой жидкости

(1)

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

2. Основные положения численного метода

В одномерном случае и = и I, Ур = рх I, У и = их и т.д., где I - направляющий вектор оси х. С учетом этого систему (1) можно записать следующим образом

+ и • У/7'к = вк . (2)

Здесь Т к, - компоненты векторов

1 = (р, и, р)т, О = (-рУи, -р"1рх, -рС32У-и)т.

Неизвестными функциями наряду с р, и, р считаются также и их пространственные производные рх, их, рх. Они определяются не дифференцированием функций р, и, р, а из уравнений

х + и •У/х' = о' - их • У/' , (3)

полученных дифференцированием уравнений (2).

Как и в работах [7, 9], система (2), (3) расщепляется на конвективную и неконвективную части. Конвективная часть имеет следующий вид

(4)

fk + u-Vfk = 0 , fk + u-Vff = 0,

имеют вид рг + и •Ур = -рУ • и,

—1

иг + и •Уи = —р Ур, р{ + и •У р = - р С§ У • и.

Здесь t - время, р - плотность, и - скорость, р давление, С5 = у]Г(р + В)/р - скорость звука, Г, В - кон- Здесь X = х - х,

а неконвективная - следующий ? = вк , х = о' - их • У/'. (5)

В алгоритме используются две разнесенные сетки: сетка 50 с узлами х1 и сетка 51 с узлами х/+1/2 = хI + Ах /2 . Сеточные аппроксимации функций р, р определяются на сетке Б0, сеточные аппроксимации и - на сетке 51.

Вычисления ведутся с шагом по времени Аtn = tn+1 - tn, п - номер временного шага. На каждом шаге по времени сначала рассчитываются уравнения конвективной части (4), так что по из-

£ ' П £ ' П

вестным значениям ТI ' , т^’ искомых сеточных функций на временном слое tп находятся промежуточные значения *. После этого рассчиты-

ваются уравнения неконвективной части (5). Здесь с использованием Т'’11, 'п , ',* находятся значения

/к,п+1 , /к,п+1 искомых сеточных функций на временном слое ^+1.

Решение уравнений конвективной части. Для расчета уравнений конвективной части (4) применяется один из вариантов полулагранжевых методов [9]. Для искомых функций, определенных на основной сетке 50, имеем

ff• = Rk(x, - и"Ut"), fk• = R( - u"U")

dx

где Rf (x) - интерполяционная функция вида Rf (x) = (1 + op X)1 ^CmXm

(6)

0<m<3

станты.

a =1 при fk;" . ff:"p < 0, a = 0 при fk;" . fk;"p > 0,

коэффициенты Ст , р вычисляются из условий ^(х,) = Г'п, ^(хШр)= Т«ирп,

(х ) = /',п дЯ,- ( ) = /',п

дх \Ли х I ’ дх \ЛШр) х ,ир

где /ир = / - 1 при и, > 0, /‘ир = / + 1 при и, < 0. При этом

полагается, что С3 = 0 при а = 1, а также в = 0 при а =

0. При а = 0 функция (х) становится полиномом

третьего порядка, в противном случае это рациональная функция.

Значения скорости в узлах, где эти компоненты не определены, рассчитываются либо с использованием соответствующих интерполяционных функций вида (7), либо простым осреднением

ип = (и,п-1/2 + и/+1/2 ) / 2 . (8)

Аналогичные (6), (7) выражения применяются и для искомых функций, определенных на сетке 51.

При использовании (6) шаг по времени Аtп ограничен условием

( л \

< тіп

/

Ах

I < |

(9)

Решение уравнений неконвективной части. Уравнения = С' неконвективной части (5) аппрок-

симируются по времени следующими неявными соотношениями рп + 1 - р*

р^ =- р*У • ип+1,

А,п

п+1 *

и^ = -(р* )1Урп+1,

рп+1 - Р*

А,п

(10)

= -р*с:;2у-и

.п+1

При такой аппроксимации схема расчета р/

п+1

и,

п+1

V'

рп+1 - р* р*С;*2А^2

+ -

V- и*

А^п

п+1

(11)

начинается с опреде-

Расчет р,п+1, и/п+1, р, ления р, из уравнения (11). Пространственная дис-

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

„ п+1

вых значений давления р/ решается итерационным

методом последовательной верхней релаксации. Затем

/1 II п+1

из второго уравнения (10) находится и/ , а

_ п+1

р/ определяется из соотношения

рп+1 = р* + (рп+1 - р* )(С5/)-2.

Уравнения х = Ох' - и х •ут'

неконвективной части (5) для искомых производных, определяемых на основной сетке 50, аппроксимируются следующими соотношениями

£ кп+1 — £ к * 'хі —'хі

Єк /■'к

і+1 — Ьі -

1 — £ к * +1 —1 хі 2Ах ’

(12)

АГп 2 Ах "

где С^ = (Т+п+1 - //+*)/А,п . Такая аппроксимация для С±1 основывается на том, что уравнение = С' неконвективной части (5) аппроксимируется следующим выражением (Т'п+1 - *)/А,п = С' .

Соотношения (12) применяются в конце алгоритма расчета неконвективной стадии, когда рассчитаны все значения . Как отмечалось выше, значения компонент вектора скорости в узлах, где эти компоненты не определены, рассчитываются либо с использованием соответствующих интерполяционных функций вида (7), либо простым осреднением (8). Аналогичные (12) выражения применяются и для искомых производных, определяемых на сетке 51.

Искусственная вязкость. Введение искусственной вязкости сводится к добавлению в правые части второго и третьего уравнений системы (1) слагаемых Ои и Ор, соответственно

V, Ур Л

и, + и У и = —— + Ои, р

р{ + и •Ур = -рС| У • и + Ор, где

1

Ои = — у я», Ор = -(Г - 1)д„ у • и, р р

(13)

Яч - счрI — СЭ Аи +

Г +1 2

(14)

рГ ’ не имеет ограничений на шаг по времени. Входящие в (10) пространственные производные приближаются центральными разностями на разнесенных сетках 50, 51.

Применив ко второму уравнению системы (10) операцию дивергенции и выразив У • ип+1 из третьего дующим образом уравнения, легко получить

Аи = тт(0, их Ах),

оу - коэффициент искусственной вязкости [10].

Таким образом в алгоритм расчета добавляется еще одна стадия - учета искусственной вязкости. Если полученные на неконвективной “невязкой” стадии значения скорости и давления обозначить и^!, и р"И+1,к (сеточный индекс / для краткости

опущен), то применяемые на стадии учета искусственной вязкости соотношения можно записать сле-

,п+1

- ип;1 + А^аи

рп+1 = рп+1 + аго; ,

ои , Ор определяются через значения сеточных

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

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

А,п = К

Ах

СНТ

тах(| ип | +С£)

(15)

где Кент - число Куранта, кСНТ < 1. Метод настоящей работы позволяет проводить расчеты и при КСНТ > 1 при выполнении ограничения (9). Однако при кСНТ > 1 может возникать чрезмерное размывание крутых волн. Численные эксперименты показали, что оптимальное значение Кскл- для расчета ударных волн находится в интервале 0,1 < КСЕЛ- < 1. При КСЕЛ- < 0,1 погрешность практически не уменьшается, а при КСЕЛ- > 1 она с ростом КСЕЛ- довольно быстро увеличивается.

3. Результаты расчетов

Исследование особенностей расчета ударных волн в переменных плотность, скорость, давление с помощью разработанной методики на основе метода С1Р-СиР выполнялось на одномерных задачах о столкновении двух одинаковых по плотности и давлению потоков жидкости, движущихся навстречу друг другу с одинаковыми (по модулю) скоростями. В результате такого столкновения в жидкости образуются две симметричные плоские ударные волны, распространяющиеся в нормальном к плоскости столкновения направлении противоположно движению потоков. Параметры задачи задаются следующим образом: в начальный момент времени t = 0 жидкость (вода с параметрами Г = 7,15, В = 3072 бар) с плотностью р = 103 кг/м3 и давлением р =1 бар движется со скоростью и = -и0 в области х > х0 и скоростью и = и0 в области х < х0, где и0 > 0. Полагается, что х0 = 30 м, и0 = 50 и 500 м/с. Расчетная область 0 < х < 60 м покрывается равномерной сеткой из 300 ячеек. Расчеты проводились с шагом по времени, выбираемым согласно (15) при Корт = 0,2.

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

При введении слагаемых с искусственной вязкостью в уравнения для скорости и давления (13) и правильном выборе коэффициента с» в численном решении на ударных волнах будут выполняться соотношения Гюгонио, которые выводятся из законов сохранения на разрывах массы, импульса и полной энергии. Если слагаемые с искусственной вязкостью добавляются только в уравнение для скорости (т.е. Ор = 0), то

при правильном выборе коэффициента с» в численном решении на ударных волнах будут выполняться соотношения, получаемые из условий сохранения на разрывах массы, импульса и энтропии. Энтропия на ударных волнах, вообще говоря, не сохраняется, но ее изменение является величиной порядка £3, где £ = (р2 -р1 )/р1 < 1, р1, р2 - плотности перед и за фронтом ударной волны соответственно. С учетом этого, для ударных волн с £3 << 1, что в жидкости может быть и

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

При относительно слабом столкновении потоков (со скоростью и0 = 50 м/с), когда изменение энтропии на возникающих ударных волнах мало (имеет порядок £3 и 3-10-5), применение искусственной вязкости не является необходимым (рис. 1). В этом случае в численном решении при корт = 0,1 и Су = 0 осцилляции на фронтах ударных волн, распространяющихся в обе стороны от разрыва, отсутствуют. Давление за ударными волнами и положение их фронтов хорошо согласуются с аналитическими решениями, полученными как в условиях сохранения на скачке полной энергии, так и в предположении постоянства энтропии (в обоих случаях здесь и далее подразумевается, что кроме этих условий на скачке выполняются также и условия сохранения массы и импульса). При и0 = 50 м/с эти аналитические решения графически совпадают.

Рис. 1 - Давление в правой половине расчетной

^-^400 п

области в момент t и > а, — 0,01 с. Ско-

/-^п=1 ’

рость сталкивающихся потоков 50 м/с. Жирная линия - графически совпадающие аналитические решения с сохранением на разрыве полной энергии и энтропии, тонкая линия с символами • - численное решение без искусственной вязкости

При увеличении скорости сталкивающихся потоков до и0 = 500 м/с (рис. 2) образуются ударные волны с изменением энтропии, которое уже нельзя считать малым (порядка £ и 0,016). При этом в аналитическом решении, полученном при условии сохранения на разрыве энтропии, скорость распространения ударных волн и давление за их фронтом оказывается заметно ниже, чем в решении с сохранением полной энергии. Это объясняется тем, что при использовании условия сохранения энтропии имеющее место на ударной волне преобразование механической энергии в тепловую не учитывается. В численном решении без введения искусственной вязкости в окрестности фронта удар-

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

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

Рис. 2 - Давление в правой половине расчетной об-

V -\800 п

ласти в момент t и ^ а, — 0,01 с. Скорость

сталкивающихся потоков 500 м/с. Жирные линии - аналитические решения с сохранением на разрыве полной энергии (линия 1) и энтропии (линия 2), тонкая линия с символами • - численное решение без искусственной вязкости

При использовании искусственной вязкости в уравнениях для скорости и давления (13) и определении Ои , Ор согласно (14) в численном решении, рассчитанном при Су = 0,8, осцилляции на фронтах ударных волн отсутствуют (рис. 3). При этом положение самих фронтов и давление за ними удовлетворительно согласуются с аналитическим решением, полученным при выполнении на разрыве закона сохранения полной энергии. При меньших значениях Су согласование с аналитическим решением ухудшается, что на рис.3 демонстрирует кривая, соответствующая Су = 0,2.

“Г

^0 2 Рис. 3 - Давление в правой половине расчетной обу -(800 п

ласти в момент t и ^ а, = 0,01 с. Скорость

сталкивающихся потоков 500 м/с. Жирная линия -аналитическое решение с сохранением на разрыве полной энергии, тонкие линии с символами - численные решения, полученные при определении Ои , Ор согласно (14) с еу — 0,2 (*) и еу — 0,8 (•)

При использовании искусственной вязкости только в уравнении для скорости (т.е. при О; =0) и с

определением Ои согласно (14) численное решение, рассчитанное при Су = 0,8, удовлетворительно согласуется с аналитическим решением, полученным при

Рис. 4 - Давление в правой половине расчетной

4 -(800 п

области в момент t и > А, = 0,01 с. Ско-

4—1п =1

рость сталкивающихся потоков 500 м/с. Жирная линия - аналитическое решение с сохранением на разрыве энтропии, тонкие линии с символами - численные решения, полученные при Ои согласно (14), Ор —0 с еу — 0,2 (х) и еу — 0,8 (•)

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

Заключение

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

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

При расчете ударных волн с существенным изменением энтропии на фронте без введения искусственной вязкости в окрестности фронта ударной

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

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

Ширина размывания фронта ударных волн при их расчете методом С1Р-СИР соответствует противо-потоковым схемам первого порядка точности.

Работа выполнена в рамках программы РАН

и при поддержке РФФИ (проект 12-01-00341-a).

Литература

1. О.С. Кочнева, Г.И. Павлов, Ж.М. Сахабутдинов, Вестник Казан. технол. ун-та, 1, 137-143 (2003).

2. А.С. Галушко, Р.Ш Абиев, Вестник Казан. технол. унта, б, 199-20З (2008).

3. А.П. Панов, Ультразвуковая очистка прецизионных деталей. Машиностроение, Москва, 1984, 88 с.

4. Р.Р. Набиев, К.А. Терещенко, Н.В. Улитин, Вестник Казан. технол. ун-та, 15, 21, 9З-9б (2012).

З. Ф.М. Гимранов, А.Н. Беляев, И.В. Флегентов, А.С. Суслов, Вестник Казан. технол. ун-та, 8, 289-291 (2012).

6. S. Kami, J. Comput. Phys., 111, 1, 31-43 (1994).

7. T. Yabe, P.Y. Wang, J. Phys. Soc. Japan, 60, 7, 210З-2108 (1991).

8. А.А. Аганин, Т.С. Гусева, Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки, 154, 4, 74-99 (2012).

9. T. Yabe, F. Xiao, T. Utsumi, J. Comput. Phys., 169, 2, ЗЗб-З9З (2001).

10. Y. Ogata, T. Yabe, Comput. Phys. Commun, 119, 2-3, 179-193 (1999).

© А. А. Аганин - д.ф.-м.н., проф., зав. лаб. Вычислительной динамки сплошной среды ИММ КазНЦ РАН, [email protected]; Т. С. Гусева - к.ф.-м.н., н.с. лаб. ВДСС ИММ КазНЦ РАН, [email protected]