Численное Решение задач взаимодействия жидкости и газа на эйлеровых сетках без явного выделения контактных границ Текст научной статьи по специальности «Физика»

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

УДК 519.6:532:533 Т. С. Гусева

ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧ ВЗАИМОДЕЙСТВИЯ ЖИДКОСТИ И ГАЗА НА ЭЙЛЕРОВЫХ СЕТКАХ БЕЗ ЯВНОГО ВЫДЕЛЕНИЯ КОНТАКТНЫХ ГРАНИЦ

Ключевые слова: контактная граница, сжимаемые среды, эйлерова сетка, неявное отслеживание контактной границы,

метод CIP-CUP.

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

Keywords: contact interface, compressible fluids, Eulerian grid, interface capturing, CIP-CUP method.

A numerical technique has been realized for computing problems of contact interaction between compressible fluids with strong deformations of the contact boundary and shock waves. It is based on the method CIP-CUP (Constrained Interpolation Profile - Combined Unified Procedure). Computations are conducted on the Eulerian mesh without explicit separation of contact boundaries. Comparison of numerical solutions of some test problems with their known solutions is presented, showing operability of the created algorithm. The performance capabilities of the realized technique has been illustrated by computing problem of impact of axially-symmetric high-speed liquid jet upon a rigid wall and a thin layer of a liquid on the rigid wall.

Введение

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

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

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

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

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

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

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

воздействии на стенку кумулятивной струйки жидкости, иллюстрирующие возможности методики.

1. Математическая модель

Рассматривается задача взаимодействия жидкости и газа при больших изменениях положения и формы контактной границы. Динамика сред без учета эффектов вязкости и теплопроводности описывается следующими уравнениями р, + и • V р = - рУ- и,

-1

и, + и • V и = -р Ур,

I К Н, (1)

р{ + и • V р = - р 01 У-и, ф, + и -Уф = 0.

Здесь

0Э = ф 0Э, I + (1 _ ф) 0Э, д ,

0з,1 =у1 Г(р + В)/р , 05,д =у1 ур/р ,

, - время, р - плотность, и - скорость, р - давление, 0э,, 0э,д - скорости звука в жидкости и газе, Г, В -константы уравнения состояния Тэта, у - показатель адиабаты, ф - функция-идентификатор среды. Предполагается, что 0 < ф < 1: в области жидкости ф = 1, в области газа ф = 0, а в малой окрестности их контактной границы функция ф непрерывна и монотонно меняется от 0 до 1 .

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

При численном решении входящего в (1) уравнения переноса для функции-идентификатора ф могут возникать ее чрезмерная диффузия и немонотонность в окрестности контактной границы. Для того чтобы избежать этого, решается уравнение р + и• УР = 0 ,

в котором Р = 1д[(1 - е) п (ф - 0,5)] [13], е - параметр, определяющий ширину области значений ф с ненулевым градиентом (ширину размывания контактной границы), полагается, что £ = 0.01. После определения Р посредством обратного преобразования находится ф.

С учетом этого систему (1) можно записать в следующем компактном виде

+ и • У1к = вк . (2)

Здесь fk, Ок - компоненты векторов 1 = (р, и, V, р, Р)т,

О = (-рУи, -р_1рх, -р_1ру, -р0з2У-и, 0)т.

В методике расчета неизвестными считаются функции fк и их производные по пространственным координатам: f кх, f ку. Они рассматриваются как независимые искомые функции и определяются не по узловым значениям функций fк, а из следующих уравнений, полученных дифференцированием уравнений (2)

fk + и • VfXk = ОХк - их -У^, fk + и ^Ук = вку - иу-У^.

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

ftk + и • = 0, fXkt+ и •VfXk = 0, fyу + и • Vf к = 0 (4)

и неконвективную часть

ftk = Gk, fk = Gkx - ux • Vfk, $ = Gky - uy • Vfk. (5)

Решение уравнений конвективной части. Для расчета уравнений конвективной части (4) применяется метод RCIP [11] из группы CIP (Constrained Interpolation Profile) полулагранжевых методов. В соответствии с методом RCIP для искомых функций на конвективной стадии имеем fk *

'ii 1 'ij \"ij

k

fk•= Rj (x, _ unДtn ).

fk • = fxij =

dR,

(x , - unДtn ( - uJAt"),

dx

dy

где Rj (x) - интерполяционная функция

Rk(x) = (1 + аюРюX + a0#0Y)-1 ^C^XrYs . (6)

(7)

0<Г+ в<3

Здесь X = х - X/, У = у - у/,

0,0 = 1 при ’У,? • fXУ,inUp] < 0, Ою = 0 при ’У,? • fXУ,inUp] > 0,

а01 =1 пРи ’у/ • у, и < с1, а01 = 0 при УП • fyk;■ и > 0

коэффициенты 0Г5, рю, Р01 вычисляются из условий

Як (х) = Ч в узлах ;*up/, ^/up, ;upjup,

УЯу (х) = У’у,п в узлах ц, /ир/, /р,

где ^ = ;- 1 при и;/ > 0, ^ = / + 1 при и;/ < 0, /^ = / - 1 при V/ > 0, /ф = / + 1 при V/ < 0. При этом полагается, что 03о = 0 при а10 = 1 и 0о3 =0 при а01 = 1, а также рю = 0 при а10 = 0 и р01 = 0 при а01 = 0. При а10 = а01

= 0 функция Як (х) становится полиномом третьего

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

Шаг по времени Д,л при расчете конвективной части ограничен условием

' ^ ^ -1 ^

ДҐ < min

ij

I un j

1 1 Дx Ду

(В)

Решение уравнений неконвективной части. Уравнения ’у = Gk неконвективной части (5) аппроксимируются по времени следующими неявными соотношениями

3n+1 _ р^ = _ діпр^ •un+1,

u ' - W = -Діп (р^ )-1Vp

n+1

і

n+1

n+1

(9)

p" ‘ - p^ =^tnpCs2V un+1.

Наряду с этим полагается Рл+1 = Р*, поскольку неконвективная часть в уравнении для Р отсутствует. При

П+1 П+1 П+1

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

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

V

pn+1 _ p^ p^Cf^n2

V. u

Mn

(lO)

Пространственная дискретизация равенства (10) определяется пространственной дискретизацией уравнений (9).

Расчет уравнений неконвективной части начинается с решения уравнения (10) методом последовательной верхней релаксации. Затем по найденным р;/ "+1 из второго уравнения (9) находится и Л+1

п+1

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

р? *' = р/ +р +- р

л/j

C

•2

s.ij

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

^ к?+1 - *к* пк пк п п п ?

1X;/ ’х;/ _ Gi+1/ Gi_1/ . к* и;+1/ и;-1/ . к* vi +1/ vi-1/

Mn

yij

yij

2Дx

*»Л /-\Л

2Дx

yij

2Дx

Діп

Gk /~*k ..n ..n ..n ..n

= ij +1 _ Gij _1 _ f j • U/j+1 _ U/j _1 _ f j • vij+1 _ vij _1

= 2Ду xij 2ДУ yij 2ДУ .

2Ду kij 2Ду

С к /хкл+1 1к * \ < «хп . „

;±1/' = ('/±1/ - ч±1/) / Д и аналогично для / ± 1.

Искусственная вязкость. Схема СГР-СИР не является консервативной. Поэтому при решении задач с интенсивными ударными волнами требуется введение искусственной вязкости. При этом к описанным выше конвективной и неконвективной стадиям расчета добавляется стадия учета искусственной вязкости. Если полученные на неконвективной “невязкой” стадии значения скорости и давления

обозначить ипн+]зк и рп+зк (сеточные индексы У, /

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

n+1 n+1

+ Д^, pn+1 = pnH+L + MnQp.

где

Qu =—Vqv, Qp = -(к- 1)qv V.u ,

P

qv = cv р |- Cs Ди + — Ди21, Ди = min(0, AV. u)

2

К = ф Г + (1- ф) Y , Су = ф Су,, + (1-ф) Суд ,

Су,/, Суд - коэффициенты искусственной вязкости жидкости и газа, Л - параметр, имеющий порядок характерного размера ячеек расчетной сетки [10].

Выбор шага по времени. В настоящей работе шаг по времени выбирается следующим образом

1

11

(11)

тах(| и? | +0?/) ^ Дx Ду, где Ксят - аналог числа Куранта, 0,1 < к0ят ^ 1.

3. Тестовые расчеты

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

Вращение кругового профиля с вырезом. Для тестирования алгоритма решения уравнений переноса рассмотрена плоская задача о переносе профиля в циркуляционном потоке (рис. 1) [12]. Изначально профиль представляет собой возмущение с единичной амплитудой нулевого значения функции ф в круговой области с вырезом. Задача решается в безразмерных переменных. Сетка равномерная 100*100 ячеек с шагом Дx = Ду = 1, шаг по време-

x

ij

невязк

ни Дt = 0,01, компоненты скорости определяются выражениями и = ы (у - ус), V = - ы (X - xc), где угловая скорость ы = 1, координаты центра вращения Xс = ус = 50.

В рамках контактного взаимодействия эту задачу можно интерпретировать как движение в поле постоянной скорости несжимаемой среды с круговым включением с вырезом, занятым другой несжимаемой средой, при отсутствии вязкости и поверхностного натяжения. Значения 0 и 1 функции ф идентифицируют области, занятые разными средами. Граница между средами находится там, где 1 < ф < 0. Согласно начальной дискретизации функция ф изменяется от 0 до 1 в области шириной в одну ячейку. Поэтому в поле изолиний ф граница возмущенной области в начальный момент t = 0 (начальная граница между средами) выглядит как четкий контур. На рис. 1 приведены также изолинии ф после одного оборота вокруг центра вращения, рассчитанные при использовании интерполяционной функции (6) с возможностью переключения между кубическим полиномом и рациональной функцией. Видно, что профиль функции ф после одного оборота почти полностью совпадает с начальным, то есть граница возмущенной области (граница между средами) остается четкой. Это достигается за счет тангенциального преобразования Р = 1д[(1 - е) п (ф - 0,5)]. Однако не для всех схем оно приводит к повышению эффективности расчета контактных границ. В частности, оно не устраняет значительных дисперсионных и диффузионных погрешностей в схемах Лакса-Вендроффа и “против потока” соответственно [13].

Рис. 1 - Профиль возмущения (изолинии функции ф) в начале вращения и после одного оборота (f1 = 2п/ш), рассчитанного методом RCIP

Взаимодействие акустической волны с контактной границей. Для тестирования алгоритма расчета неконвективной части (5) рассматривалась одномерная задача о взаимодействии акустической волны с границей раздела сред типа газ-жидкость [14]. Задача решалась в безразмерных переменных. Газ находится в области 0 < х < 0,5, жидкость - в области 0,5 < х < 1. Для газа Y = 1,4, для жидкости Г = 7,15, B = 3072,04. В начальный момент времени газ и жидкость покоятся (u = 0) при давлении p = 1, плотность газа р = 1,025-10"' 3, плотность жидкости р = 1. Акустическая волна входит в область газа через ее левую границу x = 0, где задается p = 1 + 0,1 (1 - cos ы t), ы = 2nCSg0/A, Cs,g° = 37 - невозмущенная скорость звука в газе, А = 0,225 - длина волны. Расчеты проводились при Дх = 1,2510"3, Дt = 510"6 на сетке из 800 ячеек. Как и в работе [14], вычисления выполнены без учета конвективных слагаемых.

На рис. 2 видно, что акустическая волна в ходе распространения по газу практически не размывается (штриховая кривая). В процессе взаимодействия с границей раздела между газом и жидкостью она частично отражается, а частично проходит через нее. Результаты настоящей работы хорошо согласуются с результатами работы [14].

21

1

Л \ /

/\ у ^

* 1 1 •— і

0.25 0.5 0.75 X 1

Рис. 2 - Взаимодействие акустической волны с границей контакта газа и жидкости (отмечена вертикальной пунктирной линией). Распределения давления в момент, когда волна еще не дошла до контактной границы (линия 1), и в ходе взаимодействия волны с контактной границей (линия 2). Символы • - результаты работы [14]

4. Иллюстрация применения методики

Для иллюстрации возможностей разработанной методики рассматривается задача об ударном воздействии на стенку кумулятивной струйки жидкости, возникающей на поверхности кавитационного пузырька при его схлопывании вблизи стенки. Изучение воздействия таких струек на стенку представляет значительный интерес в связи с проблемой кавитационного разрушения поверхностей. Расчеты [15] показывают, что наибольшие скорости струйки достигаются в случае, когда ее диаметр мал по сравнению с размерами пузырька. Поэтому влиянием верхней и боковых границ пузырька можно пренебречь, и ударное воздействие кумулятивной струйки можно моделировать как удар столбика жидкости по стенке [16].

Рассматривается осесимметричная задача об ортогональном ударе цилиндрической струи жидкости с плоским концом по неподвижной плоской твердой стенке или по тонкой жидкой прослойке на ней. Струя считается бесконечно длинной и окруженной неограниченным объемом газа. Эффекты вязкости, теплопроводности жидкости и газа, а также влияние поверхностного натяжения не учитываются. Для жидкости начальные значения параметров полагаются следующими: р = 103 кг/м3, Г = 7,15, В = 3047 бар, ф = 1; а для газа - следующими: р = 1 кг/м3, Y = 1,4, ф = 0. В начальный момент времени давление всюду равно 1 бар. Радиус струи Я = 25 мкм, ее скорость 500 м/с. Эти значения соответствуют кумулятивной струйке жидкости, образующейся на поверхности кавитационного пузырька в воде при комнатных условиях в том случае, когда пузырек примыкает к телу и в начале схлопывания имеет вид эллипсоида вращения с отношением полуосей ~0,84. Коэффициент искусственной вязкости

в жидкости оч,і = 0,4, а в газе Оч,д = 0,9, шаг по времени определяется из (11) при Ксет = 0,8.

Удар струи жидкости с плоским концом по твердой стенке. На рис.3 представлены результаты расчета с использованием реализованной в работе методики в области (0 < г/Й < 1,4)*(0 < ї/Я < 1) на равномерной расчетной сетке 400*300 ячеек.

Рис. 3 - Ударное воздействие струи с плоским концом на стенку через ~6 нс после начала удара. Изолинии давления: 1 - ударная волна, 2 - волна разрежения; изолинии функции-идентификатора ф: 3 - контактная граница (граница струи), пунктирная кривая - начальное положение контактной границы

В настоящей работе граница струи определяется неявно, ее положение, рассчитанное на эйлеровой сетке, показано жирной кривой 3. Эта граница образована набором близкорасположенных изолиний функции-идентификатора ф, подобно тому, как кривая ударной волны образуется близкорасположенными изолиниями давления. Следует отметить, что ударная волна 1 , возникающая в струе в момент удара по стенке, является довольно интенсивной. Скорость ее распространения составляет ~1900 м/с (для сравнения, скорость звука в невозмущенной струе ~1500 м/с), а давление и плотность за ее фронтом превышают невозмущенные значения соответственно в ~104 и ~1.26 раз. Форма струи, положение ударной волны, значение давления за ее фронтом удовлетворительно согласуются с результатами работы [16], в которой наличие окружающего газа не учитывается, сетка является адаптивно-подвижной и подстраивается под границу струи, а уравнения движения жидкости решаются с применением модификации метода С. К. Годунова второго порядка точности.

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

На рис.4 представлены результаты расчета удара струи жидкости с плоским концом по жидкой прослойке на твердой стенке в области (0 < г/Я < 1,8)*(0 < г/Я < 1) на сетке 550*300 ячеек. Как видно, данный случай характеризуется сильной деформацией

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

Рис. 4 - Ударное воздействие струи с плоским концом на тонкий слой жидкости на стенке через ~6 нс после начала удара. Обозначения кривых как на рис. 3

Заключение

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

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

Автор выражает признательность профессору А. А. Аганину за полезные обсуждения.

Работа выполнена в рамках программы РАН и при поддержке РФФИ (проект 12-01-00341-а).

Литература

1. О.С. Дмитриева, А.В. Дмитриев, А.Н. Николаев, Вестник Казан. технол. ун-та, 16, 3, 63-65 (2013).

2. В.А. Алексеев, С.В. Алексеев, Р.И. Насрыев, Р.Р. Ахме-

тов, Вестник Казан. технол. ун-та, S, 182-185 (2011).

3. L.A. Crum, J. de Physique, 40, 8, C8-285 - C8-288 (1979).

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

5. S. Kami, J. Comput. Phys., 112, 1, 31-43 (1994).

6. R. Abgrall, S. Kami, J. Comp. Phys., 169, 2, 594-б23 (2001).

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

8. T. Yabe, F. Xiao, T. Utsumi, J. Comput. Phys., 169, 2, 55б-593 (2001).

9. F.H. Harlow, A.A. Amsden, J. Comput. Phys., 3, 1, 80-93 (19б8).

10. Y. Ogata, T. Yabe, Comput. Phys. Commun., 119, 2-3, 179-

193 (1999).

11. F. Xiao,Month. Weath. Rev., 128, 4, 1165-1176 (2000).

12. S. T. Zalesak, J. Comput. Phys, 31, 3, 335-362 (1979).

13. T. Yabe, F. Xiao, Computers Math. Applic., 29, 1, 15-25 (1995).

14. M. Ida, Comput. Phys. Commun., 150, 3, 300-322 (2003).

15. А.А. Аганин, М.А. Ильгамов, В.Г. Малахов, Т.Ф. Халитова, Н.А. Хисматуллина, Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки, 153, 1, 131-146 (2011).

16. А.А. Аганин, М.А. Ильгамов, Т.Ф. Халитова, В сб. Актуальные проблемы механики сплошной среды. К 20-летию ИММ КазНЦ РАН. Т.1. Фолиант, Казань, 2011. C.134-146.

© Т. С. Гусева - к.ф.-м.н., н.с. лаб. Вычислительной динамки сплошной среды ИММ КазНЦ РАН, ts.guseva@mail.ru.