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

____________УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА

Том 156, кн. 2 Физико-математические науки

2014

УДК 519.6:532:533

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

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

Аннотация

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

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

Введение

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

55

56

А.А. АГАНИН, Т.С. ГУСЕВА

с преобразованием координат порядок точности метода существенно снижается даже при небольшом нарушении однородности сетки [4].

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

1. Методика расчета

1.1. Построение сетки. В двумерном случае soroban-сетка представляет собой набор узлов, расположенных на линиях, параллельных, например, оси Ox (рис. 1). Количество узлов на разных линиях может быть различным. На каждом временном шаге строится новая сетка узлов xn+1 = (xn+1 ,yn+1), адаптирующаяся к решению и никак не связанная со старой сеткой узлов xn = (xn, yn): линии могут перемещаться в поперечном направлении, в случае, приведенном на рис. 1, -в направлении оси Oy, узлы на каждой из линий могут перемещаться вдоль нее, количество линий и узлов может изменяться. Для нумерации линий и узлов вводятся независимые одномерные массивы, узлы нумеруются сквозным образом.

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

M(x,t)

1 + а

dx

2

+ в

df

d2x ,

x е [xs,xe].

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

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОНТАКТНОГО ВЗАИМОДЕЙСТВИЯ... 57

8 9 10

4 5

^-3

/=1

-7=1

Рис. 1. Двумерная soroban-сетка, j - номер линии, i - номер узла

Введя функцию

X

I(x,i)_j и(,М) <!.,

можно определить координаты сетки из Nn+1 узлов на временном шаге n + 1, исходя из условия, что приращение функции I на каждом отрезке будущей сетки постоянно и равно I(xE)/ (Nn+1 — 1):

^n+1 i — 1 f i 1 (xE )

Nn+1 - 1 '

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

xn+1

xi

(II

n

i+1

I0)xn + (I0 — In)xn+1 ТП ТП

Ii+1 Ii

где In _ I(x'n,tn), Io £ [Ii,Ii+1 ] • Количество узлов сетки может выбираться произвольным образом. В данном алгоритме оно задается как округленное до целого числа отношение интеграла монитор-функции по всей длине линии к заданному максимальному шагу сетки

N "+1 ^(АЧ) + 1.

\ ^^max /

Соотношение максимального и минимального шагов можно ограничение на монитор-функцию

U(x,t) _ min I у 1 + а

+в

d2f

д2 x

Um

регулировать, вводя

x £ [xs, xe],

где Umax - константа. Если U принимает везде минимальное значение (равное 1), то строится равномерная сетка с шагом Axmax. В противном случае минимальный шаг сетки будет определяться как Axmin _ Axmax/Umax.

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

U (x,y,t)

\

1 + а I ( тт-

dx

+ , f

9Уу

+ в

d2f д 2f

д2 x + d2y

2

x £ [xs,xe], y £ [ys,УЕ]•

58

А.А. АГАНИН, Т.С. ГУСЕВА

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

Хе

My(y,t) = q max M(x,y,t) + (1 — q)---------- M(x,y,t) dx, 0 ^ q ^ 1.

xs^x^xE Xe — xs J

xs

Затем для каждой новой линии yn+1 определяются новые позиции узлов с использованием значений функции Mj (x,t) = M(x,yn+1,t), получаемых линейной интерполяцией.

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

др

dt

др

dt

+ u ■Ур + u ■ Ур

р У ■ u,

ди

dt

pcS у ■ u,

+ u ■ yu

др

dt + u y

Ур

р

0.

(1)

Здесь

Cs = pCsi + (1 — p)Cs2, Csi = л/Гг(р + Bi)/p,

p - плотность, u - скорость, p - давление, Csi - скорость звука, ri,Bi - константы уравнения состояния Тэта для жидкости, ri = y, Bi =0 для газа ( y -показатель адиабаты), р - функция-идентификатор среды. Предполагается, что 0 ^ р ^ 1: в области первой среды р =1, в области второй среды р = 0, а в малой окрестности их контактной границы р непрерывна и монотонно меняется от 0 до 1 .

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

I + uyf = G, (2>

где f = (р, ^р,р), G = (—рУ ■ u, —р-1Ур, —pC2sУ ■ u, 0).

Как и в работах [1, 4-6], численное решение системы (2) разбивается на конвективную и неконвективную стадии

г»=Н

где f

f * f n

Atn

f n+1 _ f

Atn

+ u yf =0,

G,

результат расчета конвективной стадии (3).

(3)

(4)

1.4. Решение уравнений конвективной стадии. Для расчета группы уравнений переноса (3) применяется метод RCIP (Rational-Cubic Interpolation Propagation), являющийся одним из вариантов метода CIP. Решение уравнения переноса

в узле xn+1 в момент tn+1

| + uyf = 0

tn + Atn можно приближенно записать в виде

f* = fn (xn+1 — u(xn+1)Atn),

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОНТАКТНОГО ВЗАИМОДЕЙСТВИЯ... 59

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

В методе RCIP кроме узловых значений самой функции f независимыми искомыми параметрами являются также производные функции в узлах сетки, в двумерном случае - по обеим координатам fx, fy. Далее дифференцирование по x и у обозначается как дх(),ду() соответственно, а для обозначения производных, которые предполагается аппроксимировать конечными разностями, используются обозначения д()/дх,д()/ду. Неизвестные fx, fy определяются уравнениями, полученными дифференцированием (2) по х и по у соответственно:

f + U ■Vfx = dxG - dxU -Vf, f + u -Vfy = dy G - dy u ■ Vf.

Решение этих уравнений также осуществляется путем расщепления на конвективную стадию

fx - f

+ u ■ Vfx

Atn

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

f f n

fy f y

Atn

+ u ■Vfy = 0,

fx = fx - Atn (uxf* + vxfy), fy = fy - Atn (uyfx + Vyf*),

0

(5)

(6)

и неконвективную стадию

fn+1 - fx = dG fn+1 - fy = dG

Atn dx, Atn dy

(7)

В соотношениях (6) используются значения ux, uy, vx, vy, которые находятся одновременно с u(xn+1).

Одномерная интерполяционная функция, определенная на отрезке [xi,xi+1], имеет вид [6]

й-С1Р(Уд fi+1, fxi, fxi+1, X, Xi+1 xi)

(1 + abX)-1 CmXm

0<m<3

Здесь a =1 при fxifxi+1 < 0, a = 0 при fxi fxi+1 > 0; коэффициенты Cm, b вычисляются из условий равенства функций f (x) и fx(x) в узлах xi, xi+1 заданным значениям fi, f+1, fxi, fxi+1. В случае a = 0 интерполяционная функция становится полиномом третьей степени, а в случае a =1 полагается C3 = 0, так что она представляет собой отношение полинома второй степени к полиному первой степени.

Если x0 = x"n+1 - u(x"n+1)Atn - координата отправной точки, принадлежащей отрезку [xn,xn+J старой сетки, то решение одномерного уравнения переноса в точке новой сетки x"n+1 в момент tn + Atn определяется как

f *(xn+1) = RCIPfif+fif

-n X xn — xn)

xi+1, X,xi+1 xi),

где X = x0 - xn. Значение производной можно получить дифференцированием интерполяционной функции по координате x

f*(xn+1) = dxRCIPf, fn+1, x, f^i+1,X, xn+1 - xn).

60

А.А. АГАНИН, Т.С. ГУСЕВА

(xn + 1 , yn +1) a k Nm 6 N Np

<- •—j1 >(x y0) У Cm\ ' Sm c ; C S Sp

/1 -A 'i , /1+1 ~ C - - - Л- 12

S

O--------------X------------* O'-------------------X

Рис. 2. (а) Отправная точка (x0,y0) на soroban-сетке. Сплошные линии и символы ■ -новая сетка шага tn+1, штриховые линии и символы • - старая сетка шага tn , символами х и о обозначены точки, значение функции f в которых определяется интерполяцией; (б) особенности дискретизации на soroban-сетке: узел C - центральный узел шаблона. Штриховой линией показан контрольный объем, используемый для расчета давления в центральном узле

В двумерном случае расчет конвективной части осуществляется с использованием расщепления по направлениям Ox, Oy. Пусть (x°, y°) = (xn+1 — u(xrn+1)Atn,

уП+ — v (хП+ )Atn) - координаты отправной точки, в которой требуется рассчитать значение функции f и ее производных по известному их распределению в момент tn на старой сетке хп (рис. 2, a). Для расчета на двух линиях старой soroban-сетки таких, что yni ^ у0 ^ yn2, j2 = j1 + 1, отыскиваются пары узлов таких, что ХП1 ^ х0 ^ хП1+1 и хП ^ х0 ^ хП2+1. Далее, значения функции f и ее производных определяются в точках Sfas = х0) на линии j1

fs = RCIP(fn, fn+i, fni, fni+i, х0 — хП1, хП1+1 — хП1), fxS = dxRCIP(fH, fn+i fii ,f”i+i ,х0 — хП1, х”+1 — хП1) и N (хn = х0) на линии j2

fN = RCIP(fi2, fi2 + 1, fXni2, fXni2 + 1, х0 — хП2, х£ + 1 — х^),

fxN = dxRCIPf, fi2+1, fxi2, fni2+1, х° — хП2,хПг + 1 — х^) .

Для получения f(х0,у0) путем одномерной интерполяции на отрезке [ys,yN] нужны значения fys и fyN:

fyS = RCIP(fynii,fyni1 + 1, fynxii, fyxii + 1, х0 — хП1, хП1 + 1 — хП1), fyN = RCIP(fyni2,У2+1, fyxi2, fnxi2+1, х0 — х^^ — хn).

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

fyxS = dxRCIP(fy

fyxN = dxRCIPf

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

f

dy

d2G

0 n n — хП1)

1,х xi1, xi1 + 1

0 n n — хП2)

1, х xi2, xi2+1

fy

fix — Atn

^ дх

uy Q__ + uyxfx + ux fy x + vy fy x + vx Q_ + vyxfy

f n+1 _ f

yx —

(8)

yx

Atn

дyдx

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОНТАКТНОГО ВЗАИМОДЕЙСТВИЯ... 61

Финальная интерполяция на отрезке [ys,yN] определяет решение

f*(xni+1,ynj+1) = RCIPfs, fN, fys, fyN, yU - y”i,yj2 - y”i)

№

n+1

*3

j1

) = dyRCIP(fs, fN, fys, fyN, yU - yjl ,y% - y”l)

f*x(xn+17y:n+1) = RCIP(fxs, fxN, fyxs, fyxN, yU - yni,yn2 - y^), f*xy (xn+1,yn+1) = dy RCIP(fxs, fxN, fyxs, fyxN, yU - - yj1).

Расчет конвекции выполняется уже на новой сетке (отправные точки определяются для узлов новой сетки x”+1), а значения f, fx, fy, fxy в отправной точке определяются интерполяцией с использованием fn, fn, ffny на старой сетке xn. Таким образом, здесь фактически объединены перемещение сетки и перенос значений со старой сетки на новую.

1.5. Решение уравнений неконвективной стадии. Система (4) записывается в виде следующих неявных соотношений:

pn+1 - р*

Atn

un+1 - u* Atn

pn+1 - p*

At"

р* У ■ un+1,

ypn+1

р*С*2У ■ un+1.

(9)

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

pn+1 - p* ^ fVpn+1\ У ■ u*

p*Cs2Atn2 = V р * ) Atn

(10)

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

u

*+1

u

-Atn

(11)

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

Ap

р* C*s2Atn2

f yAp

У ■ u*+1 Atn

(12)

где Ap = pn+1 - p*. Система линейных алгебраических уравнений относительно узловых значений Ap, получаемая при дискретизации (12), решается методом последовательной верхней релаксации. Затем корректируется значение скорости

un+1 = u*+1 - Atn yAp (13)

р

и вычисляются давление

on+1 = p* + Ap,

(14)

а также плотность с учетом того, что согласно третьему уравнению из (9) y^un+1

= - Ap/ Atnр* C*b?,

р

n+1

* + ap р + CW.

(15)

62

А.А. АГАНИН, Т.С. ГУСЕВА

1.6. Особенности численной реализации. Построение конечно-разностных аппроксимаций производных на неструктурированной soroban-сетке имеет некоторые особенности, подробно описанные в [5, 7]. Производные по x в узле C (рис. 2, б) аппроксимируются центральными разностями на неравномерной сетке с использованием узлов C, Cp, Cm . Для построения конечно-разностных аппроксимаций производных по у в узле C требуются значения расчетных параметров в точках S и N (рис. 2, б). Для их определения применяется RCIP-интерполяция.

Узлы soroban-сетки не связаны между собой ячейками, поэтому модифицированный метод CIP-CUP применяется на неразнесенной сетке, когда давление, плотность и компоненты вектора скорости вычисляются в одном и том же узле. С целью повышения устойчивости и точности расчетов на неразнесенной сетке конечно-разностная аппроксимация некоторых слагаемых строится подобно аппроксимации на разнесенных сетках [5]. В частности, аппроксимация дивергенции скорости в узле C (рис. 2, б) для уравнения (12) строится с использованием значений компонент скорости на гранях прямоугольника, проходящих через середины интервалов между рассматриваемым узлом и соседними узлами. Аппроксимация отношения Vp/p в узле C для уравнений (11), (13) строится посредством линейной интерполяции с использованием значений Vp/p на гранях этого прямоугольника. При этом значения параметров в серединах интервалов между узлами сетки также находятся с использованием RCIP-интерполяции.

Вместо построения конечно-разностных аналогов dG/dx, dG/dy, d2G/dydx при аппроксимации уравнений (7) и второго уравнения из (8) в силу того, что влияние G уже учтено при вычислении f n+1 на неконвективной стадии, используются соотношения [5, 7]

f n+1

J X

fX +

d(fn+1 dx

f)

fn+i 7 , d(fn+1 - f)

fy = fy + dy

а затем

f n+1 J yx

fyx +

d(fn+1-

dx

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

стадия

un+1

u

n+1

невязк

+ Atnq„,

n+1

p

n+1

невязк

+ AtnQp,

где un++K и рП+вЯзк - значения, полученные на неконвективной «невязкой» стадии

(11)-(15),

Qu Vqv, Qp (к 1)qvV * u

P

qv = CvP (-Csди + ^ДИ2) , ди = min(0, Л V * u),

К = ур1 + (1 - у+2, Cv = <fiCv,1 + (1 - <fi)Cv,2,

cv,1, cv 2 - коэффициенты искусственной вязкости сред 1 и 2, Л = (ДxДy)1/2, Дx, Ду - локальные шаги сетки [5].

В задачах с ударными волнами при использовании явных схем, например схемы

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

Дx

Atn = +rt min

(\un \ + Cn)

(16)

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОНТАКТНОГО ВЗАИМОДЕЙСТВИЯ... 63

а в двумерном случае - по формуле

Atn

кскт min

1 J_

(|un| + CS) \Ax + Ay

1

(17)

где min[ ] - минимальное значение аргумента во всей расчетной области, кскт -число Куранта, кскт < 1 • Хотя условие Куранта не имеет отношения к схеме CIP-CUP, оно будет применяться в дальнейшем. В [3] показано, что при использовании метода CIP-CUP приемлемыми являются значения числа Куранта 0.1 < кскт < 1, при кскт > 1 может возникать чрезмерное размывание фронта ударной волны.

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

Разработанные алгоритмы и программы расчета протестированы на ряде задач в одномерной и двумерной постановках. Сравнивались численные решения, полученные с использованием реализованной ранее методики на основе метода CIP-CUP на фиксированных разнесенных равномерных сетках (эту методику будем называть старой) и новой методики на основе модифицированного метода CIP-CUP на неразнесенных сетках: динамически адаптивных soroban-сетках и фиксированных равномерных сетках.

2.1. Соударение движущейся и покоящейся частей жидкости. В области, занятой жидкостью с параметрами Г = 7.15, B = 3072 бар, ро = 103 кг/м3, ро = 1 бар, в начальный момент времени t = 0 имеется разрыв скорости:

\ио, х > хо, u = <

10, x < x0,

хо = 35 м, uo = —200 м/с.

На рис. 3 представлены точное решение данной задачи о распаде разрыва и результаты расчетов, полученные на фиксированных равномерных разнесенной (а-в) и неразнесенной (г-е) сетках из 100 узлов с шагом Ах = 70/100 м, а также на soroban-сетке с Axmin = Ах/2 и Axmax = 7Ах (ж-и). Построение soroban-сетки производилось при a = 3, b =1.5. Начальное распределение параметров при расчете на soroban-сетке задавалось на равномерной сетке из 200 узлов. В таком случае начальная область разрыва, к которому адаптируется сетка, имеет меньшую ширину, чем если бы начальные данные задавались на сетке из 100 узлов, что положительно сказывается на качестве решения. Расчеты проводились при cv = 0.3, кскт = 0.1 в (16) на равномерных сетках и кскт = 0.2 на soroban-сетке. Шаг по времени во всех случаях составлял примерно 3.6 • 10~5 с.

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

64

А.А. АГАНИН, Т.С. ГУСЕВА

Р/Ро\

103-

1.05- Г v6 \ 0- u/\u0\- п

Р/Ро- -0.5- \

L_ 1- 1 —1 1 Ц -1-

1

a

в

0

Рис. 3. Пространственные распределения давления, плотности и скорости в момент t = = 0.015 с в задаче о соударении покоящейся и движущейся частей жидкости, полученные с применением старой методики на равномерной разнесенной сетке (а-в), по новой методике на неразнесенных сетках: равномерной сетке (г-е) и soroban-сетке (ж-и). Точное решение - кривые без точек, численные решения - кривые с точками (значениями в узлах соответствующей сетки)

Рис. 4. Зависимость количества узлов soroban-сетки от времени в ходе расчета соударения покоящейся и движущейся частей жидкости

образом, что шаг сетки оказывается здесь в 2 раза меньше шага используемой равномерной сетки. В результате фронты ударных волн становятся более крутыми. Кроме того, решение в пределах размытого фронта ударной волны, распространяющейся влево, и в области прохождения ударных волн становится монотонным, не возникает провала плотности в окрестности точки соударения. Количество узлов soroban-сетки на протяжении всего расчета не превышает 53 (рис. 4), в то время как равномерная сетка, обеспечивающая аналогичное разрешение фронта ударных волн, состояла бы из 200 узлов.

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОНТАКТНОГО ВЗАИМОДЕЙСТВИЯ... 65

Р>Ро

104-

0

P/Po 104 -

1.2-

P/Po

1.1-

1-

rx

0-т

0-

u/u*

1 x/x0 2

1

1 x/x 0 2

0

3

и

1

0

0

Рис. 5. Пространственные распределения давления, плотности и скорости в момент t = = 0.0092 с в задаче о распаде разрыва в жидкости, полученные с применением старой методики на равномерной разнесенной сетке (а-в), по новой методике на неразнесенных сетках: равномерной сетке (г-е) и soroban-сетке (ж-и), и* = 264 м/с. Точное решение -кривые без точек

2.2. Распад разрыва в жидкости. В покоящейся жидкости (Г = 7.15, B = = 3072) в начальный момент времени t = 0 имеется разрыв параметров:

{14089ро, х ^ xo, I 1.216ро, x ^ хо,

Р = \

Ро, х < хо, [ро, х < хо,

Ро = 1 бар, ро = 103 кг/м3, хо = 20 м.

На рис. 5 представлено точное решение этой задачи и результаты численных расчетов на таких же равномерных сетках и при тех же параметрах построения soroban-сетки, что и в предыдущей задаче. Начальные поля плотности, давления и скорости при расчете на soroban-сетке также задавались на равномерной сетке из 200 узлов. Расчеты проводились при cv = 0.3, &crt = 0.1 в соотношении (16) на равномерных сетках и cv = 0.5, &crt = 0.2 на soroban-сетке. Шаг по времени во всех случаях составлял примерно 2.3 • 10~5 c.

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

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

66

А.А. АГАНИН, Т.С. ГУСЕВА

Рис. 6. Зависимость количества узлов soroban-сетки от времени в ходе расчета распада разрыва в жидкости

Рис. 7. Пространственные распределения давления, плотности и скорости при t = 1 с в задаче о распаде разрыва на границе газ-жидкость, полученные с применением старой методики на равномерной разнесенной сетке (а-в), по новой методике на неразнесенных сетках: равномерной сетке (г-е) и soroban-сетке (ж-и), и* = 540 м/с. Точное решение -кривые без точек

полученное на soroban-сетке решение ближе к точному, а используемое при этом количество узлов на протяжении всего расчета остается меньше 100 (рис. 6).

2.3. Распад разрыва на границе газ —жидкость. В начальный момент времени область x > xo занята жидкостью (р = 1) с параметрами Г = 7.15, B = 3072 бар, p = 14089p0, р = 1.216 р0, а область x < x0 - газом (р = 0) с параметрами y = 1.4, p = po, р = 10~3 • ро . Обе среды неподвижны, po = 1 бар, р0 = 103 кг/м3, x0 = 1200 м.

На рис. 7 представлено точное решение задачи о распаде разрыва на границе газ - жидкость и результаты расчетов, полученные в области длиной 5000 м с применением фиксированных равномерных разнесенной (а-в) и неразнесенной (г-е) сеток из 400 узлов с шагом Ax = 5000/400 м, а также на soroban-сетке с Axmin =

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОНТАКТНОГО ВЗАИМОДЕЙСТВИЯ... 67

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

= Дж/2 и Дxmax = 7Дx (ж-и). Начальное распределение параметров при расчете на soroban-сетке задавалось на равномерной сетке из 500 узлов. Расчеты проводились при cv,i = 0,cvg = 0.7. Для равномерной сетки &crt = 0.4 в соотношении (16), для soroban-сетки &crt = 0.8. Шаг по времени в обоих случаях составлял около 2 • 10~3 с.

В начальный момент времени влево по газу начинает распространяться ударная волна, вправо по жидкости - волна разрежения. На рис. 7, б, д, з приведена часть расчетной области x/xo < 2 для того, чтобы более детально показать ударную волну в газе и следующую за ней контактную границу. В отличие от предыдущих задач, здесь новая методика на равномерной сетке (рис. 7, г-е) дает в целом удовлетворительное решение. Только положение фронта ударной волны, как можно видеть на рис. 7, д, оказывается менее точным по сравнению с остальными двумя решениями. При использовании soroban-сетки за счет задания начального разрыва на более мелкой сетке и своевременного сгущения сетки все особенности решения описываются более точно. Как видно по вставкам на рис. 7, а, г, ж, более детально изображающим решение в окрестности ударной волны в газе и заднего фронта волны разрежения в жидкости, при использовании soroban-сетки давление жидкости за контактной границей (при t =1 с ее положение x/xo ~ 0.5) оказывается более близким к точному решению: задний фронт волны разрежения отделяется от контактной границы, в отличие от решений, полученных с использованием равномерных сеток. В решении, полученном на soroban-сетке, положения ударной волны и контактной границы значительно лучше согласуются с точным решением, чем в обоих случаях равномерных сеток, фронт ударной волны становится заметно круче. Вместе с тем в результате разрежения в областях перед и за фронтом ударной волны в газе, между контактным разрывом и задним фронтом волны разрежения, а также между передним и задним фронтами волны разрежения количество узлов soroban-сетки на протяжении всего расчета остается меньше 200 (рис. 8).

2.4. Цилиндрический взрыв в неоднородной среде. В начальный момент времени в цилиндрической области взрыва (x — 0.5)2 + (у — 0.5)2 < 0.08 находится газ с у = 1.2, у = 0, p = 3, р =1 ,а вне ее - газ с у = 1.8, у = 1, p = 0.1, р = 0.1. Оба газа неподвижны. В начальный момент времени на границе области возникают радиально расходящаяся ударная волна и радиально сходящаяся к центру области взрыва волна разрежения. На рис. 9 приведены профили плотности перед началом (1) ив процессе (2) фокусировки волны разрежения

68

А.А. АГАНИН, Т.С. ГУСЕВА

Рис. 9. Задача о цилиндрическом взрыве в неоднородной среде. Радиальные профили плотности в сечении у = 0.5 в моменты t = 0.03 (кривая 1) и 0.075 (кривая 2)

в центре области взрыва. Профили рассчитаны с применением старой методики на фиксированной разнесенной сетке 401 х 401 (Ax = Ay = 1/400), kcRT = 0.4 в (17).

На рис. 10 представлены результаты расчетов этой задачи на момент t = 0.075 с применением стационарных разнесенных и неразнесенных равномерных сеток, а также soroban-сетки. Использовались равномерные сетки с 101 х 101 узлами c Ax = Ay = 0.01 и с 201 х 201 узлами c Ax = Ay = 0.005, soroban-сетка с Axmin = = 0.01/2, Axmax = 10 • 0.01 (a = 0.3, b = 0.52, q = 1). Расчеты проводились при kcRT = 0.4 в соотношении (17) на равномерных сетках 101 х 101 и &crt = 0.8 на равномерных сетках 201 х 201 и на soroban-сетке. Шаг по времени во всех случаях составлял примерно 5.5 • 10~4. Применялась искусственная вязкость с cv = 0.6.

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

При использовании soroban-сетки начальное распределение параметров в расчетной области задавалось на равномерной сетке с 201 х 201 узлами. Применение густой начальной сетки и, соответственно, более точное задание начальной границы разрыва, к которой адаптируется soroban-сетка, заметно улучшает круговую симметрию численного решения. Изолинии плотности на рис. 10, з получены интерполяцией решения с soroban-сетки на равномерную сетку 201 х 201. На рис. 10, ж приведен вид soroban-сетки в момент t = 0.075. Сетка довольно резко сгущается в окрестности фронта ударной волны, контактной границы и в областях с локальными экстремумами (в окрестности минимальных значений плотности сходящейся волны разрежения). Динамическая адаптация soroban-сетки позволяет получить более качественное, чем в обоих случаях равномерных сеток 101 х 101, решение: круговая симметрия улучшается, фронт расходящейся ударной волны становится более крутым, а минимумы плотности в области сходящейся волны разрежения - более выраженными. Кроме того, при использовании soroban-сетки

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОНТАКТНОГО ВЗАИМОДЕЙСТВИЯ... 69

0.8

0.6

0.4

0.2

0.2 0.4 4 0.6 0.8

1

0.8 -

0.6 -

>4

0.4 -0.2 -0

t/ГС

0 0.2 0.4 X 0.6 0.8 1

Рис. 10. Задача о цилиндрическом взрыве в неоднородной среде. Верхний ряд - результаты, полученные по старой методике на разнесенных равномерных сетках, средний ряд -по новой методике на неразнесенных равномерных сетках, нижний ряд - по новой методике на soroban-сетке; (а-д, з) - изолинии плотности в момент t = 0.075, внутренняя жирная кривая - начальная граница области взрыва, штриховая кривая - текущее положение контактной границы (р = 0.5); (в, е, и) - радиальные профили плотности в сечении у = 0.5 в момент t = 0.075: точками показано решение, полученное на сетках 101 х 101, символами + - на сетках 201 х 201 и на soroban-сетке. Точки и символы + на профилях плотности соответствуют узлам сеток. Вид soroban-сетки в момент t = 0.075 (ж)

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

В данном случае применения soroban-сетки число ее узлов не превышает количества узлов равномерной сетки 101 х 101 на протяжении всего расчета (рис. 11) за счет сильного разрежения сетки в областях, где решение изменяется мало. Время расчета в рамках новой методики оказывается на soroban-сетке в 1.8 раз меньше, чем на равномерной сетке 201 х 201, а полученное решение оказывается столь же

70

А.А. АГАНИН, Т.С. ГУСЕВА

200x200-1

100x100-

0-1----------------------1---

0 t 0.07

Рис. 11. Зависимость количества узлов soroban-сетки от времени в ходе расчета цилиндрического взрыва в неоднородной среде

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

Заключение

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

Работа выполнена при поддержке РФФИ (проект № 14-01-97004р_поволжье_а).

Summary

A.A. Aganin, T.S. Guseva. Numerical Simulation of the Dynamics of Non-Uniform Compressible Media Based on the CIP-CUP Method on Dynamically Adaptive Soroban Grids.

We compare the efficiency of the CIP-CUP method on a stationary uniform staggered grid and the modifications of this method on a non-staggered dynamically adaptive Soroban grid and a non-staggered stationary uniform grid in their use for computing problems of wave propagation in uniform (liquid) and non-uniform (with the presence of gas-gas and gas-liquid contact boundaries) media. We show that as compared with the initial version of the CIP-CUP method on a stationary uniform staggered grid, its modification on a Soroban grid allows one to conduct computations with similar or even higher accuracy and to apply significantly smaller number of grid vertices. The modification of the CIP-CUP method on a stationary uniform non-staggered grid is least efficient. Moreover, it can lead to non-monotonic solutions.

Keywords: shock waves, non-uniform medium, CIP-CUP method, adaptive Soroban grids.

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОНТАКТНОГО ВЗАИМОДЕЙСТВИЯ... 71

Литература

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

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

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

4. Yabe T., Mizoe H., Takizawa K., Moriki H., Im H.-N., Ogata Y. Higher-Order Schemes with CIP Method and Adaptive Soroban Grid towards Mesh-Free Scheme // J. Comput. Phys. - 2004. - V. 194, No 1. - 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. - V. 40, No 1. - P. 167-183.

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

7. Аганин А.А., Гусева Т.С. Методика расчета волн в жидкости и газе методом CIP-CUP с применением динамически-адаптивных Soroban-сеток // Вестн. Башкир. ун-та. - 2014. - Т. 19, № 2. - С. 368-380.

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

24.03.14

Аганин Александр Алексеевич - доктор физико-математических наук, профессор, заведующий лабораторией, Институт механики и машиностроения КазНЦ РАН, г. Казань, Россия.

E-mail: aganin@kfti.knc.ru

Гусева Татьяна Сергеевна - кандидат физико-математических наук, старший научный сотрудник, Институт механики и машиностроения КазНЦ РАН, г. Казань, Россия.

E-mail: ts.guseva@mail.ru