Методика расчета волн в жидкости и газе методом CIP-CUP с применением динамически-адаптивных Soroban-сеток Текст научной статьи по специальности «Физика»

УДК 519.6:532:533

раздел ФИЗИКА

МЕТОДИКА РАСЧЕТА ВОЛН В ЖИДКОСТИ И ГАЗЕ МЕТОДОМ CIP-CUP С ПРИМЕНЕНИЕМ ДИНАМИЧЕСКИ-АДАПТИВНЫХ SOROBAN-СЕТОК

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

Институт механики и машиностроения Казанского научного центра РАН Россия, Республика Татарстан, 420111 г. Казань, ул. Лобачевского, 2/31.

Тел./факс: + 7 (843) 236 52 89.

E-mail: immkazan@mail.ru

Разработаны алгоритмы и программы расчета двумерных волн в неоднородных средах без явного выделения контактных границ методом CIP-CUP с применением динамически-адаптивных Sornban-сеток. Проведено их тестирование и апробация на одномерных и двумерных задачах. Полученные численные решения сравнивались с аналитическими решениями и численными решениями других авторов. Проводилось сравнение численных решений, полученных на Soroban-сетках и стационарных равномерных декартовых сетках. Показано, что в ряде случаев Soroban-сетки позволяют при значительно меньшем числе узлов получить численное решение, более близкое к точному не только по величине, но и по характеру изменения.

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

Введение

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

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

В своих первоначальных вариантах метод CIP-CUP был ориентирован на использование стационарных декартовых расчетных сеток. Один из таких вариантов с применением разнесенных сеток был реализован и авторами [2, 3]. Однако в расчетах с ударными волнами и быстро изменяющимися контактными границами более предпочтительны адаптивные сетки, сгущающиеся в окрестности больших градиентов решения и разреженные в областях, где градиенты решения малы. С учетом этого в работе [4] был предложен новый вариант метода CIP-CUP, в котором дискретизация расчетной области проводится с помощью Soroban-сетки (от японского "soroban" - счеты). Это неструктур иро-ванная сетка, узлы которой не связаны между собой никакими элементами типа ячеек. С одной сто-

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

В настоящей работе на основе изложенной в [4] методики разработаны алгоритмы и программы расчета двумерных волн в неоднородных средах без явного выделения контактных границ с применением динамически-адаптивных 8огоЪап-сеток. Проведено их тестирование на одномерных и двумерных задачах, а также сравнение с численными решениями, полученными на равномерных стационарных декартовых сетках.

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

В двумерном случае 8огоЪап-сетка определяется рядом параллельных оси x прямых линий и расположенным на них множеством узлов (рис. 1). В ходе расчета линии могут перемещаться в направлении оси у и их количество может изменяться. Узлы на каждой из линий могут перемещаться в направлении x, их количество может изменяться и быть различным для разных линий. На каждом временном шаге, фактически, строится новая сетка, никак не связанная со старой. Таким образом обеспечивается гибкость в адаптации сетки. Для нумерации линий и узлов вводятся независимые одно-

* автор, ответственный за переписку

мерные массивы. Нумерация узлов осуществляется сквозным образом. Номер узла непрерывно увеличивается от линии к линии и уникален для всей сетки.

Рис. 1. Вид двумерной 8огоЬап-сетки, jL - номер линии, iG - номер узла.

С целью определения областей, где необходимо сгущение сетки, вводится монитор-функция, которая в одномерном случае определяется как линейная комбинация длины дуги и кривизны кривой, являющаяся зависимостью какого-либо параметра решения f (например, плотности) от пространственной координаты x,

M (x,t) = min

1 +

+ ß

dx2

, M m

x £ [xs , xE ] ,

Здесь Мтах - константа, позволяющая регулировать соотношение максимального и минимального шагов будущей сетки, xS, xE - координаты начала и конца разбиваемого отрезка. В областях с возрастающими градиентами решения функция M(x,t) также возрастает. Введя функцию

л

I (л, t) = J M (л,

можно определить координаты адаптивной сетки на временном шаге п+1 из Ln+1 узлов, исходя из условия, что приращение функции I на каждом отрезке новой сетки - постоянная величина, равная !(хЕ)/(£п+1-1),

хп+1 _ г-1| . 1 (ХЕ )

Х _ IV1-1

Здесь I"1 - обратная к I функция, значения которой могут быть получены линейной интерполяцией

лГ1 =

(I»! - I о) л» + (10 - I» ) xh

тП тП

Ii+1 - Ii

где Ii = I(xnk, tn), 10 e[/f, /гв+1]. Количество узлов сетки может выбираться произвольным образом. В данном алгоритме оно определяется как округленное до целого числа отношение интеграла монитор -

функции по всей длине линии к заданному максимальному шагу сетки

Ln+1 = int

' «Хв)Л Ar

V^^ max j

+1.

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

Лхтт Лхтах/Мтах.

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

M (л, у, t) = 1 +

V

2 Г +

2

+ ß

Г о о \

df + df v Эл2 ду2 ,

При построении сетки сначала определяются позиции линий. Для этого используется одномерная монитор -функция

My(y,t) = q max M(x,y,t) +

+ (1 - q)-

1

л -л

E

IM (л, у, г)йл, 0 < q < 1

S х3

Затем для каждой новой линии jL определяются новые позиции узлов с использованием распределений М^ (х, t) _ М (х, у[ jL], t), получаемых линейной интерполяцией.

1.2. Уравнения динамики сжимаемой среды В случае двух контактирующих сред их динамика без учета эффектов вязкости и теплопроводности описывается уравнениями

pt + u • Vp = -pV • u,

Здесь

ut + u • Vu = -p Vp , pt + u •Vp = -pC^V • u , Ф( + u • Vф = 0.

CS = ФСХ1 + (1 - Ф)СХ 2 , CSi = Vri (P + Bi VP

(1а) (1б) (1в)

(1г)

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

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

Все уравнения системы (1) имеют вид

+

в

S

2

+

л

S

f dt

+ u-V/ = G .

(2)

Здесь / = р, и, р, ф и, соответственно, О = - р Уи, -р-1Ур, - рС/Уи, 0.

Как и в работах [1, 4, 5], уравнения (2) расщепляются на конвективную и неконвективную (для р, и, р) части

/ -/" At"

+ u-V/ = 0,

/n+1 - /* At"

= G.

(3)

(4)

1.4. Решение уравнений конвективной части

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

+ и-V/ = 0 в узле с координатами х"+1

/ = /'

в момент t n

,n+1 \

имеет вид

"+1

\udt

или приближенно

/ * = /" (

n l1xn+1 - u(x n+1)Atn),

где аргумент представляет собой отправную точку - положение трассера, который за время At n прибудет в узел новой сетки xn+1. Если положение отправной точки не совпадает с узлами сетки, значение переносимой величины в этой точке интерполируется по известным значениям / n на старой сетке xn. Так как узлы старой и новой сеток также могут не совпадать, значение u(xn+1) в узле новой сетки определяется интерполированием с применением un на старой сетке xn.

В данной работе используется один из вариантов CIP-метода - RCIP (Rational-Cubic Interpolation Propagation) [7], в котором применяется интерполяция третьего порядка в областях, где решение гладкое, и интерполяция более низкого порядка в областях с большими градиентами решения. Как будет показано ниже, кроме значений / n в методе RCIP необходимы значения производной функции в узлах сетки, в двумерном случае по обеим координатам /х , /у. Здесь они являются независимыми искомыми параметрами и определяются уравнениями, полученными простым дифференцированием (2) по х и по у, соответственно

+ u - V/x =----V/ ,

dt дх дх

-JL + u-V/y =-

dt Jy dy

dG du

dy

V/ .

(Обозначения Э()/д х, Э()/ д у в настоящей работе используются для конечно-разностных выражений.) Эти уравнения также расщепляются на конвективную часть

/х* - /

At"

+ u-V/x = 0,

(3)

fy /y

Atn

которая дополняется соотношениями

+ u-V/y = 0,

fx = /* -Atn (u"/* + v

/*),

y x J y

f -С* A " I n -С* I n -C

/y = /y -At (Un/x + Vn/1

(5)

' y

y J x

y Jy '

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

/ "+1 - f Jx Jx

Atn /"+1 - f

yy

dG

dx

dG

(6)

Ахп ду

Расчет конвективной части осуществляется с использованием расщепления по направлениям х, у и использованием одномерной интерполяции. Интерполяционная функция, определенная на отрезке [х; , хг +1], имеет вид [7]

К-С1р(Л, /г+1, /хг, /хг+1, Х, хг+1 - Х1 ) =

= (1 + авХ )-1 X СтХт .

0<т<3

Здесь а = 1 при /"•/"+! < 0, а = 0 при

/п ' /п+1 - 0 , коэффициенты Ст , в вычисляются из условий равенства значений интерполяционной функции и ее производной в узлах хг и хг +1 заданным значениям / , /¡+1, /хг, /хг+1. При этом полагается, что С3 = 0 при а = 1, а также в = 0 при а = 0. Так что при а = 0 функция ЯС1Р становится полиномом третьего порядка, в противном случае это рациональная функция - отношение полинома второго порядка к полиному первого порядка.

Значение функции на отрезке [хг , хг +1] определяется как

/ (х) = ЯС1Р(/г, /+ , /х{, /хг+1, X, Хг+1 - ), где X = х - х, а значение производной можно получить дифференцированием интерполяционной функции по координате х (здесь обозначается как дх ( ))

/х (х) = Э х ЯС1Р( /

, /г+1, /хг, /хг+1, Х, хг+1 х1) .

В двумерном случае конвекция рассчитывается следующим образом. Пусть (х0, у0)=(х п+:-и(х п+1)М п, у п+'-у(х п+1)Мп) - координаты отправной точки, в которой требуется рассчитать значение функции /и ее производных по известному их распределению в момент Iп на старой сетке х п (рис. 2). Для этого сначала определяются две линии старой 8огоЪап-сетки, такие что у[ у1] < у0 <

n

*

у[7' 1+1 ]. Затем на линиях j1 и /2=/1+1 отыскиваются пары узлов, для которых соответственно выполнены условия х[/1] < х0 < х[г1+1] и х[г2] < х0 < х[г2+1]. Далее значения функции f и ее производных определяются в точках Я на линии /1 (верхние индексы п и * опускаются)

Рис. 2. Отправная точка (х ,у ) на 8огоЬап-сетке. Сплошные линии и символы ■ - новая сетка шага t п+1, штриховые линии и символы • - старая сетка шага tп, символы х и о - точки, значение функции / в которых определяется интерполяцией.

/8 _ ЯС1Р(/й, /й+1, /хй, /хц+1, х° - хй, хп+1 - хй),

/хЯ _ д х КС1Р( fi1, /й+1, fxi1, /хЛ+1, х° - х!'1, хг'1+1 - хг'1)

и N на линии j2

fN _ КС1Р(/г2 , /,2+1, fxi2 , fxi2+1, х0 - хг'2, хг'2+1 - хг'2 ) , fxN _ дх КС1РС/г2, /,2+1, fxi2, /хг'2+1,х - хг'2,хг'2+1 - хг'2) •

Чтобы выполнить одномерную интерполяцию на отрезке [у8, удт] в направлении у, необходимы значения /у8 и Получить их можно, применяя ЯС1Р-интерполяцию

/уЯ _ КС1Р(/уЦ, /уг'1+1, /ухг'1, /ухг'1+1,х - хг'1, хг'1+1 - хг'1) , fyN _ КС1Р(/у12, /уг2+1, /ух12, /ух12+1,х0 - хг'2,х12+1 - х12) .

При этом требуется введение дополнительных неизвестных - смешанных производных /ху = /ух, для которых

/ух8 дх КС1Р( /уй , /уЛ+1, /ухП , /ухП+1, х0 - хг'1, хг'1+1 - хг'1)

/yxN _ д х КС1Р(/у,2, /у,2+1, /ух,2, /ух,2+1, х0 - хг'2, хг'2+1 - хг'2)

Для них также требуется пересчет после конвективной стадии, подобный (5), (6),

fyx fyx At

n dfx + n f * + n f * + Uy Uyx fx + Ux Jyx +

df *

+ vn f * + vn f + vn f *

y yx x yx y

V

fn+1 _ f

J yx J y.

dy

Atn

д 2G

dydx

, (7)

После интерполяции в направлении х выполняется интерполяция по у

/ (х0, у 0) _ ЯС1Р( / , ^ , /у8 , ^ ,

у0 - у/1, У/2 - У/1)

/у (х0, у0) _ Э у ЯС1Р(/ , ^ , /у8 , ^,

у0 - у/1, у/2 - у/1)

/х (х , у ) _ К.ОР(/хЯ , /xN , /ухЯ , /yxN,

у0 - у/1, у/2 - у/1)

/ху (х , у ) _ д у К-С1Р(/хЯ , /xN, /ух8 , /yxN,

у0 - у/1, у/2 - у/1)

В итоге решение имеет вид

/*(х-1, у-1) _/п(х0, у0),

/*(хп+1, уп+1) _ /уп(х0, у0), /*(хп+1, уп+1) _ /хп (х0, у0),

/ху(хп+1, уп+1) _ /пу(х0, уУ

Расчет конвекции выполняется уже на новой сетке (отправные точки определяются для узлов новой сетки х п+:), а значения /, /х, /у, /ху в отправной точке определяются интерполяцией с использованием

/п, /п, /п, /ху на старой сетке х п. Таким образом,

здесь, фактически, объединены перемещение сетки и перенос значений со старой сетки на новую.

1.5. Решение уравнений неконвективной части

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

n+1 *

Р _Р = _p*V • u n+1, Atn

u n+1 _ u* Atn

Vp

n+1

p

pn+1 n*

P _ P = _p*C*2V • u n+1. Atn "

(9а) (9б) (9в)

Применив к (96) операцию дивергенции и вы-

разив V • un+1 из (9в), легко получить

pn+1 _ p* p*C*2Atn2

= V^

Vp p

n+1 Л

V • u*

' Atn

(10)

Определяя поле скорости на стадии предиктора по известному полю давления

и*+1 _ и* -Л^ ^р Р

можно преобразовать (10) к виду

Ap p*C*2Atn

- = V •

rVAp N

*

Р J

V • u

t"

*+1

(11)

(12)

где Ар _ рп+1 - р*. Система линейн^1х алгебраических уравнений относительно узловых значений Ар

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

УДр

п+1

*+1

-Дгп

Р

(13)

и вычисляются давление и плотность с учетом того,

что согласно (9 в) V - и п+1 = - Ар/(Агпр* С*2),

рп+1 = р* + Др ,

п+1 * Др

Рп+1=Р +-£• с?

(14)

(15)

1.6. Решение уравнений неконвективной части

Первая и вторая производные по х аппроксимируются следующими конечно-разностными соотношениями (рис. 3):

Эх,

Э2 /

хСр хс

( хС хСт )

/Ср /С

+

хСр х'

+ (хСр хС )

С

/С - /Ст

хс X,

Ст у

Эх2

/Ср /С /с - /с

Ст

у хСр хс

хс хс

Рис. 3. Особенности дискретизации на 8огоЪап-сетке.

Символы • - узлы сетки, символы х - точки, значение функции / в которых определяется интерполяцией. Узел С - центральный узел шаблона. Штриховой линией показан контрольный объем, используемый для расчета давления в центральном узле.

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

/S = К-С1Р(/Sm , /Sp , /xSm , /xSp , хС - ХSm , ХSp - ХSm ) /N = КС1Р(/Nm , , /xNm , /х^ , хС — XNm ,х ^ — XNm )

В результате конечно-разностные аппроксимации производных по у имеют вид

ду _ Э2 /

ду2

1

УN - Уs

2

( ус - Уs +

УN - ус

/с - Ь

+(УN - ус)

ус - Уs

/N - /С /С - Ь

УN - Уs

УN - ус Ус - Уs

В этом случае при дискретизации уравнения (12) для аппроксимации производных д2Лр/ду2 кроме неизвестных Лрт , Лр^ , Лрш , Лр^ потребуются неизвестные производные Лрхт , Лрх^р , Лрх^т , Лрхлр. Чтобы избежать связанных с этим сложностей, вводится следующий итерационный процесс. Пусть на итерации т

(т)

Ар(т) = Ар(т-1) + 5р(т) = £ 5р

(т)

к=0

Тогда (12) можно переписать в виде

¿р

(т)

рС*2Дгп

-V-

Уёр

(т)

р

Др

(т-1)

р*С*2Дгп2

+У -

УДр

(т-1)

р*

V-и

*+1

,(16)

Дхп

где конечно-разностная аппроксимация слагаемых с уже известными Лр(т-1) строится с помощью ЯС1Р-интерполяции, поскольку

Архт-1 =ЭАр(т-1) / дх, а для слагаемых с 5р(т) ис-

пользуется линейная интерполяция

^ =

т (т) = 1

XSp - х.4

¿т =

((Xsp - Xs « + (Xs - Xsm ^> ) ,

1

(х щ XN )3рьш +

+ (^ - хХт )3

(т)

Второе слагаемое в (16) аппроксимируется как

1 др(т -дрСт) ^

V-

г V¿p(т) ^

*

р

хСр хс

рССр

хСр хс

1 зр^ -рт

\

рССт

хС - х

Ст

1

УN - Уs

г¿р№ - ¿рт -^

PРCN У N - уС

__1_ ¿рт -рт

рCs ус - уs

где

* * [ Xс +

Рсст = Р I , УС

РССр = Р

^ хС + хСр ^

-ТТ^, Ус

С

1

2

хСр - х

С

хт х

1

+

£

P*cs = р*[xc Ус + ys

2

О* _ о*[ х уС + УN pCN _ р I хС, 2

которые определяются также с помощью ЯС1Р-интерполяции.

Из-за отсутствия связи между узлами посредством организации их в ячейки в структуре 8огоЬап-сетки метод С1Р-СИР приходится реализо-вывать на неразнесенной сетке, когда давление, плотность и компоненты вектора скорости определяются в одном и том же узле. При этом конечно-разностная аппроксимация акустической части решаемых уравнений дополняется рядом особенностей, необходимых для повышения устойчивости и точности расчетов на неразнесенной сетке [5].

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

( /

ди

*+1

dx xCp xCm

*+1

xC + xCp

\ \

, Ус

v

/

dv

*+1

дУ Ум _ ys

*+11 xc + xCm

_ и +1 -C-—, Ус

v l 2 %

Wxc ,Ус±Ук Л Л

*+1[ ^ yc + ys _ v I xc,

2

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

Э/ p* dx

Э/ p* dy

(xCp _ xC ) (pC _ p*Cm ) +

(xcp _ xCm ) p(*Cm (xC _ xCm ) . * * . + (xC _ xCm ) (pCp _ pC )

(xCp _ xCm ) p**Cp (xCp _ xC )

(Ум _ Ус) (pc _ Ps ) + (yN _ ys) pCs ( Ус _ ys) , (ус _ ys) (pn _ pc)

(уN - уя) pcN(УN - Ус), Значения параметров в точках, не совпадающих с узлами сетки, также определяются с использованием ЯОР-интерполяции. Аналогичным образом

представляется У Ар / р* в уравнении (13).

При аппроксимации уравнений (6), (5) можно избежать построения конечно-разностных аналогов для дО/_дх, дО/ду, д2О/дудх, используя тот факт, что

Ст+1

уже учтено при вычислении / на неконвективной стадии. Из (4) следует, что

ЭО 1 Э(/п+1 - / *)

dx Atn

dx

n+1

3G_ 1 Э( fn+1 _ f •)

dy Atn Также из (6) следует, что

dy

dGy

1 Э( fyn+1 _ f*)

дх Дtn Эх С учетом этого после вычисления /п+х на неконвективной стадии применяются соотношения

/хп+1 _ /* + Э(/п+1 - /^

fr1 = f* +

dx

d(fn+1 _ f *) dy

а затем

-n+1 _ ,* , d(fyn+1 _ f*)

fn+1 = f +-

J yx Jyx

dx

1.7. Искусственная вязкость

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

Эи _ Ур — + и Уи _—^ + Ои, дt р

^ + и-Ур _-рСх2У-и + йр,

где

Ои _-1 У , йр _-(к - 1)^У- и , Р р

-С8ди + ^^ди2^, ди_тш(0,ХУ-и),

К _ фГ + (1 - ф)Г , еу _ фсу,1 + (1 - ф)с„,2, Су,1, 2 - коэффициенты искусственной вязкости сред 1 и 2, X = (Лх Лу)1/2, Лх, Лу - локальные шаги сетки.

Таким образом, в алгоритм расчета добавляется еще одна стадия

Qv = CvP

.n+1

n+1

pn+1 = pn+1 У Уневязк

+ AtnQu, + AtnQp,

п+1 п+1

где и невязк и рневязк - значения, полученные на

*

неконвективной "невязкой" стадии (11)-(15), а Ои, й*р определяются через значения сеточных функций, рассчитанных на конвективной стадии.

1.8. Выбор шага по времени

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

2

1

1

u

невязк

из известного условия устойчивости Куранта. В частности, в случае равномерной эйлеровой сетки полагают в одномерном случае

Дхп = кскт тш

Дх

(I ип I+Сп)

и в двумерном случае

Дхп = кскт т1п

1

/

(I ип I +Сп)

11

— + —

Дх Ду

-1

(17)

,(18)

где тт[ ] - минимальное значение во всей расчетной области, кскт - число Куранта, кскт < 1. Метод настоящей работы позволяет проводить расчеты и при ксдт > 1 [1, 5]. Однако при ксдт > 1 может возникать чрезмерное размывание крутых волн. Численные эксперименты [3] показали, что оптимальное значение ксдт для расчета ударных волн находится в интервале 0.1 < кСКТ < 1.

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

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

2.1. Соударение движущейся и покоящейся частей жидкости

В начальный момент в жидкости с плотностью р0 = 103 кг/м3 и давлением р0 =1 бар имеется разрыв скорости: и

Параметры жидкости: Г = 7.15, В = 3072 бар. Полагается х0 = 20 м, и0 = -500 м/с. Расчетной областью является отрезок [0, 40 м].

На рис. 4 представлены точное решение данной задачи о распаде разрыва и результаты расчетов, полученные на равномерной сетке с шагом Дх=40/300 м (а-в) и одномерной 8огоЪап-сетке с Дхтт=Дх/2 и Дхтах=2Дх (г-е). Параметры для построения 8огоЪап-сетки принимались следующими а= 3, в = 0.5. Расчеты проводились при су = 0.4, ксдт = 0.1 в (17) на равномерной сетке и ксдт = 0.2 на 8огоЪап-сетке. Шаг по времени в обоих случаях составлял 5 10-6 с. Видно, что при использовании равномерной сетки в окрестности фронта ударной волны, распространяющейся по неподвижной жидкости влево, возникают осцилляции. При использовании же 8огоЪап-сетки осцилляций нет, фронт ударных волн становится более крутым. Кроме того, уменьшается величина провала плотности в окрестности точки соударения х/х0= 1. При этом количество узлов 8огоЪап-сетки на протяжении всего расчета остается меньше 300 (рис. 5), в то время как равномерная сетка, обеспечивающая аналогичное разрешение фронта ударных волн, состояла бы из 600 узлов.

2.2. Распад разрыва в однородной жидкости

В начальный момент времени г = 0 в однородной покоящейся (и = 0) жидкости (Г= 7.15, В = 3072 бар) в плоскости х0 = 20 м имеется разрыв давления. В области х > х0: р = 14088 бар, р = 1216 кг/м3, в области х < х0: р0 = 1 бар, р0 = 1000 кг/м3.

< 0 при х > х0 и и = 0 при х < х0.

щейся и движущейся частей жидкости, полученные на равномерной сетке (а-в) и 8огоЪап-сетке (г-е). Сплошные линии без

символов - точное решение, с символами • - численные решения.

300

N

0

0 t, мс 7

Рис. 5. Изменение количества узлов Soroban-сетки при расчете соударения покоящейся и движущейся частей жидкости.

На рис. 6 представлено точное решение этой задачи и результаты расчетов, полученные с применением равномерной сетки из 300 ячеек с Ax=40/300 м (а-в) и Soroban-сетки с Axmin= Ах/2 и Axmax= 2Ax (а = 3, в = 0.5) (г-е), cv = 0.4, kCrt = 0.1 в (17) для равномерной сетки и kCrt = 0.2 для Soroban-сетки. Шаг по времени в обоих случаях составлял 4-10-6 с. Как и в предыдущей задаче, при использовании равномерной сетки на фронте ударной волны, распространяющейся влево, возникают осцилляции. При применении Soroban-сетки, количество узлов которой на протяжении расчета остается меньше 300 (рис. 7), осцилляции на ударной волне не возникают, а ее фронт становится более крутым. Кроме того, лучше разрешаются волна разрежения и контактная граница.

2.3. Распад разрыва на границе газ-жидкость.

В начальный момент времени газ (у = 1.4) и жидкость (Г = 7.15, B = 3072 бар) покоятся (и = 0), а давление на контактной границе х0 = 360 м терпит

разрыв. В области жидкости (х > х0) р = 14088 бар, р = 1216 кг/м3, ф =1, в области газа (х < х0) р = 1 бар, р = 1 кг/м3, ф = 0.

300 и

N

о

О t, мс 5.6

Рис. 7. Изменение количества узлов Soroban-сетки при расчете распада разрыва в жидкости.

На рис. 8 представлено точное решение этой задачи и результаты расчетов, полученные в области длиной 1700 м с применением равномерной сетки из 300 ячеек с Ax= 1700/300 м (а-в) и Soroban-сетки с Axmin= Лх/2 и Axmax= 2Лх (а = 3, в = 0.5) (г-е). Расчеты проводились при cv,¡ = 0, cvg = 0.7. Для равномерной сетки kCrt = 0.45 в (17), для Soroban-сетки kcrt = 0.7. Шаг по времени в обоих случаях составлял около 8 10-4 с. При использовании Soroban-сетки волна разрежения, ударная волна в газе и контактная граница описываются более точно. Количество узлов Soroban-сетки лишь в небольшом временном интервале немного превышает 300 (рис. 9).

б

п

v0 Z, 1 Л/Л0 2 I Л/Л0

Рис. 6. Пространственные распределения давления, плотности и скорости в момент t = 5.6 мс в задаче о распаде разрыва в однородной жидкости, полученные на равномерной сетке (а-в) и Soroban-сетке (г-е), и* = 264 м/с. Сплошные линии без символов - точное решение, с символами • - численные решения.

10001 100 1

р/р0! 1

4 х/х0 0 1000 ^

4 х/х

д

т

0 1 2 3 4 х/х0 0 1 2 3 4 х/х0 0 1 2 3 4 х/х0

Рис. 8. Пространственные распределения давления, плотности и скорости при г = 0.35 с в задаче о распаде разрыва на границе газ-жидкость, полученные на равномерной сетке (а-в) и ЗогоЪап-сетке (г-е), и* =540 м/с. Сплошные кривые без символов - аналитическое решение, с символами • - численное решение.

300-N

2.4. Цилиндрический взрыв в однородной среде

В начальный момент времени г=0 в покоящемся (и=0) газе (у =1.4) в цилиндрической области

0 С 0.3

Рис. 9. Изменение количества узлов ЗогоЪап-сетки при расчете распада разрыва на границе газ-жидкость.

х - 0.5)2 + (у - 0.5)2 < 0.08 давление р=3, плотность р =1, а вне этой области р=0.1, р =0.1. При г > 0 формируются ударная волна, радиально расходящаяся от границы области, и волна разрежения, сходящаяся к ее центру.

0 0.2 0.4 0.6 0.8 Рис. 10. Задача о цилиндрическом взрыве в однородной среде. Верхний ряд - равномерная декартова сетка, нижний ряд - ЗогоЪап-сетка; (а, г) - изолинии плотности (внутренняя жирная кривая - начальная граница области взрыва); (б, д) -радиальные профили плотности в трех сечениях (графически совпадают): у = 0.5, х = 0.5 и х = у,

(в, е) - результаты работы [5].

б

1

Рис. 11. БогоЬап-сетки настоящей работы (а) и работы [5] (б) в расчетах цилиндрического взрыва в однородной среде, изме-

нение количества узлов сетки (в): жирная кривая - настоящая работа, тонкая кривая - работа [5].

На рис. 10 представлены результаты расчетов этой задачи с применением равномерной декартовой сетки из 101x101 узлов с Лх=Лу=0.01 и SoroЬaп-сетки с Лхтц~ Лх/2, Лхтах~ 8Лх (а = 1.1, р = 0.7, q = 1), а также результаты Takizawa и др. [5]. Шаг по времени оставался постоянным, равным 510-4. Коэффициент искусственной вязкости су = 0.55. При использовании Soroban-сетки начальное распределение параметров в расчетной области задавалось на равномерной декартовой сетке из 201x201 узлов. Это практически не оказывает влияния на время расчета, но начальная граница разрыва, к которой адаптируется Soroban-сетка, аппроксимируется с лучшим разрешением, что заметно улучшает круговую симметрию численного решения. Для решения на Soroban-сетке изолинии и профили плотности в сечениях х = 0.5 и х = у получены интерполяцией на равномерную декартову сетку 201x201. Как видно на рис. 10г, при применении Soroban-сетки изолинии плотности становятся более близкими к круговым. Вместе с тем, полученные и на равномерной декартовой, и на Soroban-сетке радиальные профили плотности в трех сечениях у = 0.5, х = 0.5 и х = у (рис. 10б, д) графически сливаются в одну линию. Символы + в профилях плотности соответствуют узлам сеток в сечении у=0.5, поэтому по их положению можно судить об отличиях распределения в этом сечении узлов равномерной сетки и SoroЬaп-сетки. При применении Soroban-сетки фронт расходящийся ударной волны становится более крутым (рис. 10д), а минимальные значения плотности в области сходящейся волны разрежения уменьшаются (рис. 10д). Расчеты показывают, что все эти особенности возникают и в численном решении на равномерной декартовой сетке при ее измельчении. Как видно, согласование представленных на рис. 10 профилей плотности настоящей работы и работы [5], рассчитанных на равномерной сетке и на близких по структуре Soroban-сетках (рис. 11), можно считать удовлетворительным.

На рис. 11а приведена Soroban-сетка настоящей работы, соответствующая тому же моменту, что и рис. 10д, а на рис. 11б - Soroban-сетка работы [5], соответствующая рис. 10е. Видно, что сетка довольно резко сгущается в областях с большими измене-

ниями решения (в окрестности фронта ударной волны) и областях с локальными экстремумами (в окрестности минимальных значений плотности сходящейся волны разрежения), что и позволяет получить более качественное, чем на равномерной сетке, решение. При этом число узлов Soroban-сетки не превышает количества узлов равномерной сетки (101x101) на протяжении всего расчета (рис. 11 в).

2.5. Цилиндрический взрыв в неоднородной среде

В начальный момент времени в цилиндрической области д/(х - 0.5)2 + (у - 0.5)2 < 0.08 находится газ с у=1.2 (ф=0), а вне ее газ с 7=1.8 (ф=1). Оба газа неподвижны (м=0). Внутри области р=3, р=1, снаружи р=0.1, р=0.1. Коэффициент искусственной вязкости Су = 0.6.

На рис. 14 представлены результаты расчетов данной задачи на равномерной декартовой сетке из 101x101 узлов с Лх = Лу = 0.01 и Soroban-сетке с ЛхтЫ~ Лх/2, Лхтах~ 8Лх (а= 7.1, р = 0.3, q = 1) в момент г = 0.075. На рис. 15 приведена Soroban-сетка в тот же момент времени, а также кривая изменения количества сеточных узлов в ходе расчета. Шаг по времени выбирался фиксированным и равным Лí=5■10"4. При использовании Soroban-сетки начальное распределение параметров в расчетной области и положение границы раздела задавалось на равномерной декартовой сетке из 201x201 узлов. Применение Soroban-сетки приводит здесь к тем же, что и в предыдущей задаче, изменениям в численном решении по сравнению с решением, полученным на равномерной декартовой сетке. А именно, изолинии плотности становятся более близкими к круговым (рис. 14в), фронт ударной волны - более крутым (рис. 14г), минимум плотности в области сходящейся волны разрежения уменьшается (рис. 14г). Кроме того, при использовании Soroban-сетки за фронтом ударной волны появляется небольшой участок со значительно меньшим градиентом плотности (рис. 14г). Расчеты показывают, что все перечисленные особенности возникают и в численном решении на равномерной декартовой сетке при ее измельчении.

Рис. 14. Задача о цилиндрическом взрыве в неоднородной среде. Верхний ряд - равномерная декартова сетка, нижний ряд -

БогоЪап-сетка. Слева - изолинии плотности (внутренняя жирная кривая - начальное положение границы раздела, штриховая кривая - текущее положение контактной границы), справа - радиальные профили плотности в трех сечениях: у = 0.5, х = 0.5 (сплошная кривая) и х = у (штриховая кривая), символы + соответствуют узлам сетки в сечении у = 0.5.

Рис. 15. БогоЬап-сетка в расчете цилиндрического взрыва в неоднородной среде (а) и изменение количества узлов сетки (б).

Заключение

Разработаны алгоритмы и программы расчета двумерных волн в неоднородных средах в переменных плотность, скорость, давление без явного выделения контактных границ методом С1Р-СИР с применением динамически-адаптивных ЗогоЪап-сеток. Для тестирования численные решения ряда одномерных задач сравнивались с аналитическими решениями, а двумерной задачи - с численными решениями других авторов. При тестировании использовались также численные решения, рассчитанные по разработанным алгоритмам и программам с применением равномерных декартовых се-

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

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

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

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

ЛИТЕРАТУРА

1. Yabe T., Wang P. Y. Unified numerical procedure for compressible and incompressible fluid // J. Phys. Soc. Japan. 1991. Vol. 60. No. 7. P. 2105-2108.

2. Аганин А. А., Гусева Т. С. Численное моделирование контактного взаимодействия сжимаемых сред на эйлеровых сетках // Ученые записки Казанского университета. Сер. Физ.-матем. науки. 2012. Т. 154. Книга 4. С. 74-99.

3. Аганин А.А., Гусева Т. С. Расчет контактного взаимодействия сжимаемых сред без явного выделения межфазных границ // Вестник Башкирского Университета. 2013. Т. 18, №3. С. 646-661.

4. T. Yabe, H. Mizoe, K. Takizawa, H. Moriki, H.-N. Im, Y. Ogata Higher-Order Schemes with CIP Method and Adaptive Soroban Grid towards Mesh-Free Scheme // J. Comput. Phys. 2004. Vol. 194. P. 57-77.

5. 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. P. 167-183.

6. Staniforth A., Cote J. Semi-Lagrangian integration schemes for atmospheric models - A review // Mon. Weather Rev. 1991. Vol. 119. p. 2206-2223.

7. T. Yabe, F. Xiao, T. Utsumi The constrained interpolation profile method for multiphase analysis // J. Comput. Phys. 2001. Vol. 169. No. 2. P. 556-593.

8. Ogata Y., Yabe T. Shock capturing with improved numerical viscosity in primitive Euler representation // Comput. Phys. Commun. 1999. Vol. 119. No. 2-3. P. 179-193.

Поступила в редакцию 07.05.2014 г.

A TECHNIQUE OF COMPUTATION OF WAVES IN GAS AND LIQUID BASED ON THE CIP-CUP METHOD ON DYNAMICALLY ADAPTIVE SOROBAN GRIDS

© A. A. Aganin, T. S. Guseva*

Institute of Mechanics and Engineering Kazan Science Center of Russian Academy of Sciences 2/31 Lobachevsky st., 420111 Kazan, Republic of Tatarstan, Russia.

Phone: +7 (843) 236 52 89.

E-mail: ts.guseva@mail.ru

Wave processes in non-uniform media with large deformations of the contact boundaries are found, for example, in apparatus of heat and mass exchange in liquids with cavitation bubbles, in devices of homogenization of liquids by cavitation, in ultrasound surface cleaning, etc. Such problems can most naturally be computed using an approach in which all the discontinuities are calculated with their smearing over a short region. One variant of such an approach based on the CIP-CUP method (by Yabe et al) on the stationary Cartesian grids has recently been realized by the authors. In this work, we apply the adaptive unstructured Soroban grids (proposed by Yabe et al), which, in the 2D case, are some variable sets of movable vertices in some variable numbers of movable parallel straight lines. A single interpolation procedure (based on the CIP method) of computing the convective terms and the values of flow-parameters at a new grid through those at an old one simplifies the use of variable grid. Numerical algorithms and computer codes of calculating two-dimensional waves in nonuniform media without explicit separation of contact interfaces by the CIP-CUP method on the dynamically adaptive Soroban grids have been developed. These codes have been tested using one- and two-dimensional problems. The calculated numerical solutions have been compared with analytical solutions and numerical results by other authors. Comparison of the numerical solutions computed on the Soroban grids and the stationary uniform Cartesian grids showed that in some cases the Soroban grids with a significantly less grid point number allow one to get such a numerical solution that is closer to the exact one not only in its magnitude but also in its behavior.

Keywords: shock waves, contact interface, shock and interface capturing, CIP-CUP method, adaptive Soroban grids.

Published in Russian. Do not hesitate to contact us at bulletin_bsu@mail.ru if you need translation of the article.

REFERENCES

1. Yabe T., Wang P. Y. J. Phys. Soc. Japan. 1991. Vol. 60. No. 7. Pp. 2105-2108.

2. Aganin A. A., Guseva T. S. Uchenye zapiski Kazanskogo universiteta. Ser. Fiz.-matem. nauki. 2012. Vol. 154. Kniga 4. Pp. 74-99.

3. Aganin A.A., Guseva T. S. Vestnik Bashkirskogo Universiteta. 2013. Vol. 18, No. 3. Pp. 646-661.

4. T. Yabe, H. Mizoe, K. Takizawa, H. Moriki, H.-N. Im, Y. J. Com-put. Phys. 2004. Vol. 194. Pp. 57-77.

5. Takizawa K., Yabe T., Tsugawa Y., Tezduyar T. E., Mi-zoe H. Comput. Mech. 2007. Vol. 40. Pp. 167-183.

6. Staniforth A., Cote J. Mon. Weather Rev. 1991. Vol. 119. p. 2206-2223.

7. T. Yabe, F. Xiao, T. J. Comput. Phys. 2001. Vol. 169. No. 2. Pp. 556-593.

8. Ogata Y., Yabe T. Comput. Phys. Commun. 1999. Vol. 119. No. 2-3. Pp. 179-193.

Received 07.05.2014.