CC BY
46
раздел ФИЗИКА
УДК 519.6:532:533
РАСЧЕТ КОНТАКТНОГО ВЗАИМОДЕИСТВИЯ СЖИМАЕМЫХ СРЕД БЕЗ ЯВНОГО ВЫДЕЛЕНИЯ МЕЖФАЗНЫХ ГРАНИЦ
© А. А. Аганин, Т. С. Гусева*
Институт механики и машиностроения Казанского научного центра РАН Россия, Республика Татарстан, 420111 г. Казань, ул. Лобачевского, 2/31.
Тел./факс: +7 (843) 236 52 89.
E-mail: [email protected]
Разработана методика численного решения задач контактного взаимодействия сжимаемых сред типа газ-жидкость. Методика позволяет рассчитывать на эйлеровых сетках сильные деформации контактной границы и ударные волны в контактирующих средах. Положение контактных границ определяется неявно с помощью функции-идентификатора среды, описываемой уравнением переноса. В основе методики - известный метод CIP-CUP (Constrained Interpolation Profile-Combined Unified Procedure), в котором уравнения динамики сред, записанные в переменных плотность, скорость, давление, расщепляются на конвективную и неконвективную части. Сравнение численных решений ряда тестовых задач с их аналитическими решениями или результатами расчетов других авторов показывает работоспособность созданной методики, а ее возможности иллюстрируются на задаче об ударе осе-симметричной высокоскоростной струи жидкости по жесткой стенке.
Ключевые слова: контактное взаимодействие сжимаемых сред, неявное отслеживание контактной границы, метод CIP-CUP.
Введение
Наиболее сложными для численного исследования контактного взаимодействия сжимаемых сред (газа, жидкости) являются задачи, в которых границы контакта могут возникать, исчезать или претерпевать значительные деформации, а акустические импедансы (произведение плотности на скорость звука) взаимодействующих сред сильно различаются. В качестве примера можно привести задачу о схлопывании кавитационного пузырька вблизи твердой стенки, которая возникает в рамках исследования механизмов кавитационного разрушения. Ранее кавитационное разрушение связывали с разнообразными негативными последствиями, например, разрушением лопастей гребных винтов водных двигателей [1]. В последнее время это явление приобрело прикладное значение, в частности, в медицине для дробления камней в почках [2]. Как показывают физические эксперименты [3-6], для схлопывания пузырька у стенки характерен сценарий с формированием направленной к стенке высокоскоростной струйки жидкости, которая пересекает полость пузырька и бьет либо по противоположной части поверхности пузырька (если пузырек не примыкает к стенке), либо непосредственно по стенке (в случае примыкания). Этот сценарий характеризуется как наличием интенсивных ударных волн, так и большими перемещениями и деформациями (даже с нарушением топологии) контактных границ.
Некоторые обзоры используемых в настоящее время методов численного решения задач динамики сплошных сред при наличии в них контактных границ можно найти, например, в [7-11].
Считается, что для задач с деформирующимися контактными границами наиболее эффективными являются лагранжевы методы. В этом случае координатные поверхности всегда совпадают с контактными границами и естественным образом деформируются в соответствии с деформацией частиц среды. Но при больших изменениях формы частиц среды лагранжева сетка сильно искажается, что приводит к потере точности и росту погрешностей. В таких случаях приходится перестраивать расчетную сетку и выполнять интерполяцию параметров ячеек со старой сетки на новую. Эти проблемы существенно уменьшаются при использовании произвольных лагранжево-эйлеровых методов. В этих методах лишь некоторые координатные поверхности (или их части) совпадают с существенными границами задачи (ударными волнами, контактными границами, стенками и т.д.). Движение остальных координатных поверхностей можно задавать произвольным образом и строить привязанные к существенным границам задачи адаптивные сетки со сгущением в областях с большими градиентами решения, с плавным изменением характерных размеров соседних ячеек, с сохранением определенной формы ячеек и т.д. Однако если существенные границы задачи, под которые подстраивается лагранжево-эйлерова сетка, сильно деформируются, возникают, исчезают или меняют свою связность, то и в этом случае приходится выполнять перестройку сетки и проводить соответствующую интерполяцию численного решения на новую сетку. Как и для лагранжевых методов это усложняет алгоритмы расчета и снижает их экономичность пропорционально частоте применения процедур перестройки сетки и интерполяции решения.
* автор, ответственный за переписку
В эйлеровых методах используются стационарные (декартовы или криволинейные) сетки. При этом подвижные контактные границы на эйлеровых сетках либо выделяются явно как совокупность поверхностных ячеек или узлов, либо отслеживаются неявно с помощью некой пространственной функции-идентификатора (например, контактная граница совпадает с линией определенного уровня этой функции или находится в области ее максимального градиента), значения которой сохраняются вдоль траекторий частиц среды. При явном выделении контактной границы в процессе расчета требуется дополнительно определять скорость перемещения узлов, составляющих эту границу. При неявном отслеживании контактной границы распределение функции-идентификатора определяется уравнением переноса с локальной скоростью среды. Методы с явным выделением контактной границы, в которых она всегда остается четкой поверхностью раздела сред, характеризуются относительно высокой эффективностью, но при обобщении на многомерный случай или при изменении топологии границ их алгоритмы значительно усложняются. При неявном отслеживании контактная граница может быть также определена как четкая поверхность раздела сред, как, например, в методах Volume Of Fluid [12] и Level Set [13]. Но соответствующие алгоритмы также усложняются при изменении топологии и увеличении пространственной размерности задачи. В то же время в рамках неявного отслеживания контактной границы существует альтернативный способ ее представления в виде тонкого переходного слоя, в котором характеристики среды (некой «смеси» контактирующих сред) плавно переходят от параметров одной среды к параметрам другой. Ширина слоя с измельчением сетки устремляется к нулю. В этом случае в качестве функции-идентификатора может выступать либо физическая характеристика среды - показатель адиабаты, постоянные двучленного уравнения состояния, например, модели «stiffened gas», либо массовая или объемная доля среды в «смеси». При этом в областях, занятых разными средами функция-идентификатор имеет постоянные различающиеся значения. В начальный момент времени функция-идентификатор может иметь близкий к ступенчатому профиль на границе контакта. В ходе расчетов вследствие численной диффузии этот профиль постепенно размывается, и точное положение контактной границы становится неопределенным. Подобная ситуация возникает при расчете ударной волны методами без явного выделения разрывов. В случае необходимости (например, для визуализации) граница контакта определяется либо как область больших градиентов функции-идентификатора, либо как изоповерхность ее среднего значения. Обусловленная численной диффузией ширина переходной области может быть ограничена при использовании определенных методик
(см., например, [14]). Методы с размыванием контактной границы являются с точки зрения описания ее положения менее точными, чем методы с четким определением ее положения, но в то же время они гораздо более просты в применении. Они легко обобщаются на многомерный случай и без дополнительных сложностей позволяют описывать контактные границы сложной топологии, включая их фрагментацию, коалесценцию, возникновение и исчезновение.
Таким образом, методы с неявным отслеживанием и размыванием контактной границы выглядят весьма привлекательными. Однако следует иметь в виду, что использование такого подхода к описанию движения контактной границы в задачах с ударными волнами, где оптимальными являются консервативные схемы сквозного счета типа схемы С. К. Годунова, оказывается проблематичным [15]. Характерная для этих схем численная диффузия, весьма желательная при расчете ударных волн в однородных средах, приводит к размыванию профиля плотности на контактной границе. В сочетании с размытым профилем функции-идентификатора на контактной границе возникает переходная зона с неопределенной средой, давление которой уже нельзя выражать через консервативные переменные (массу, импульс и полную энергию), применяя для этого одно из уравнений состояния контактирующих сред или какое-либо уравнение состояния «смеси». Такие попытки приводят к возникновению нефизичных ос-цилляций давления, которые не связаны с порядком точности схемы, поскольку наблюдаются даже при использовании схем первого порядка точности [16]. Избежать нефизичных осцилляций давления в рамках подхода с неявным отслеживанием и размыванием контактной границы можно, воспользовавшись неконсервативными схемами сквозного счета, основанными на применении уравнения эволюции давления [15, 17, 18].
Однако указанное выше сочетание способов описания движения контактной границы и расчета динамики взаимодействующих сжимаемых сред не может эффективно применяться в задачах взаимодействия газа и несжимаемой жидкости, а также в задачах, где характерное число Маха претерпевает значительные изменения. Примером является задача схлопывания кавитационного пузырька в жидкости, в которой низкоскоростная начальная стадия схлопывания является маломаховой, а в финальной высокоскоростной стадии число Маха становится выше единицы. Применение для таких задач явных методов расчета с условием устойчивости Куранта, например, метода С. К. Годунова, приводит к жестким ограничениям временного шага на маломаховой стадии. Одна из первых попыток обобщения на сжимаемые среды методов расчета несжимаемых сред, использующих уравнение для давления (или для коррекции давления), таких, как MAC [19], SMAC [20], SIMPLE [21], предпринята в методе
ICE [22]. Однако этот метод, как и его модификация PISO [23], основаны на консервативном подходе, при котором возникают трудности при расчете течений с контактными границами [14]. В относительно недавно появившемся методе CIP-CUP (Constrained Interpolation Profile-Combined Unified Procedure) [14] используются уравнения движения сжимаемой среды в неконсервативных переменных (плотность, скорость, давление), которые расщепляются на конвективную и неконвективную части. Для расчета конвективной части применяется метод CIP - один из вариантов полулагранжевых методов, характеризующийся монотонностью и малой диффузией. Для расчета неконвективной (акустической) части, учитывающей эффекты сжимаемости среды, применяется метод UP (Unified Procedure) - обобщение методов расчета несжимаемых сред, основанных на уравнении для давления. Отдельный расчет конвективных слагаемых приводит к упрощению уравнения для давления, используемого на неконвективной стадии в методе UP. В совокупности с неявным отслеживанием контактной границы метод CIP-CUP может применяться для расчета течений с широким диапазоном проявления сжимаемости сред и с сильно деформирующимися контактными границами со значительным перепадом акустического импеданса. Несмотря на то, что метод CIP-CUP не является консервативным, в сочетании с искусственной вязкостью [24] он демонстрирует хорошие возможности при расчете многофазных течений с ударными волнами. Поэтому сам метод CIP-CUP и его модификации широко применяются многими исследователями [25-32].
В настоящей работе реализована методика расчета задач контактного взаимодействия сжимаемых сред, основанная на неявном отслеживании контактной границы на эйлеровой сетке и уравнениях динамики сред в неконсервативной форме с применением для их решения метода CIP-CUP.
1. Математическая модель и методика расчета
1.1. Математическая модель
Рассматривается контактное взаимодействие между жидкостью и газом при больших изменениях положения и формы контактной границы. Динамика взаимодействующих сред без учета эффектов вязкости и теплопроводности описывается следующими уравнениями
pt + u Vp = - pV u,
ut + u Vu = -p-1Vp, t (1) pt + u -V p = -pC_2 V u,
jt + u -Vj = 0.
Здесь
Cs = jC, l + (1-j) Cs,g, CSJ =Vr(p + B)/p ,
Cs,g =VWp ,
/ - время, р - плотность, и - скорость, р - давление, С5,1, - скорости звука в жидкости и газе, Г, В - константы уравнения состояния Тэта, у - показатель адиабаты, ф - функция-идентификатор среды. Предполагается, что 0 < ф < 1, так что в области жидкости ф = 1, в области газа ф = 0, а в малой окрестности их контактной границы ф непрерывна и монотонно меняется от 0 до 1 (т.е. контактная граница находится там, где 0 < ф < 1).
1.2. Основные положения методики расчета В методике расчета вместо уравнения для ф применяется уравнение
^ + и-У^ = 0 ,
в котором ^ = гя[(1 -е)р(ф-0.5)], е - параметр,
определяющий ширину области значений ф с ненулевым градиентом. Введение функции ^ позволяет значительно сократить численную диффузию и немонотонность функции ф в области больших градиентов в окрестности контактной границы. Это особенно полезно при изменении связности контактных границ, в частности, когда границы делятся, пересекаются или сливаются.
В двумерной постановке с учетом того, что и = (и, у), Ур = (рх, ру), где х, у - оси декартовых или цилиндрических координат, систему (1) можно записать следующим образом
/к + и-У/ = Ок. (2)
Здесь /к, Ок - компоненты векторов
f = (р, и, у, р, О = (- рУ-и, - р-1рх , - р-1ру ,-рС/У-и, 0)т. В методике расчета наряду с р, и, у, р и ^ неизвестными функциями считаются также и их производные по пространственным координатам: рх, их, у х, рх, Рх, ру, иу, уу, ру, Ру. Они рассматриваются как независимые искомые функции и определяются не по узловым значениям р, и, у, р и ¥, а из уравнений, полученных дифференцированием уравнений (2) /к + и-У/к=ок - и х -У/к,
/к + и-У¡к = ок - и у -У/к (3)
Некоторые преимущества дополнения набора искомых функций р, и, у, р, ^ системы (2) их пространственными производными будут указаны ниже.
Как и в работе [14], система (2), (3) расщепляется на конвективную и неконвективную части. Конвективная часть имеет следующий вид
¡к + и-у/к=0, ¡к + и-у¡к=0,
¡к + и-У/к = 0, (4)
а неконвективная - следующий
/к = ок, /к = ох - их -у/к, ¡к=оку - и у -у/к. (5) В алгоритме расчета используются три разнесенные сетки, состоящие из прямоугольных ячеек с
координатными линиями, параллельными осям x и у: основная сетка S0 с узлами x, = (x-, у,), i, j - номера узлов вдоль осей x и у соответственно, сетка
51 с узлами (xi+1/2 , у,), xt+1/2 = x, + Ax /2 и сетка
52 с узлами (x-, у,-у v+1/2 = у v +Ду/2 . Сеточные аппроксимации функций p, р, F определяются на сетке S0, сеточные аппроксимации составляющих u и v скорости u - на сетках S1 и S2 соответственно.
Вычисления ведутся c шагом по времени At" = - t", " - номер временного шага. На каждом шаге по времени сначала рассчитываются уравнения конвективной части (4), так что по известным .. , f .. , f .. , искомых сеточных
ij ' " x, ij ' " у, ij J
функций на временном слое t" находятся промежуточные значения j, /Д*, /Д*. После этого рассчитываются уравнения неконвективной части (5). Здесь с использованием j, , и /Ц', /%
/*к,* _ .гк ,n+1 гк ,n+1 гк ,n+1 _
находятся значения / • , /x,j , искомых сеточных функций на временном слое t"+1. Такой порядок расчета, когда неконвективная стадия является заключительной, желателен в случае несжимаемой жидкости, когда необходимо обеспечить выполнение условия нулевой дивергенции скорости [28]. Для сжимаемых сред можно сначала рассчитать неконвективную часть (5), а затем -конвективную (4), как это описано, например, в работе [17].
1.3. Решение уравнений конвективной части
Для расчета уравнений конвективной части (4) применяется метод RCIP (Rational-Cubic Interpolation Propagation) [33] из группы CIP (Constrained Interpolation Profile) полулагранжевых методов. В соответствии с методом RCIP для искомых функций, определенных на основной сетке S0, имеем
дКк
" (Х„ - u"
fk• = Jj(x,-u-Dtn) fS, =it(1. -At")
f^ =—L (x - u" At"),
J yi Эу ^ . . '
(6)
где u" = (u", v"), u" = (u"-1/2, j + u-1/2, j) / 2, v" = (v" -1/2 + v"j+1/2) / 2 , Rjk (x) - интерполяционная функция
Rk (x) = (1 + 0СюP10X + 0^/)-1 S C„XTs. (7)
Здесь X = x - x,, Y = у - у,,
«10=1 при fkj • fkip. < 0, «10=0 при fkj • fk:p, > 0,
«01=1 при fkL • f. < 0,
(8)
ao1 = 0 при /kj • /j > 0, коэффициенты Д0, b01, Crs вычисляются из условий
Rk (x, )=fk ■", Rk (I*., )=fup",, К (x/, LP )=f:
Rk (x ) = f:
ij \ iup,jup J J i
k," .jup ?
/'k
iup,jup '
ЭRk , , , ЭRk
" ' 1 rk,"
(x. ) = fXkL ,1L (j ) = fx
( xlp )= fx
k," iup,j
■k," i, jup
ЭRk
(x.. ) = fk'" ,—(x. . )= f
\ . 1 J y■ij' Эу iup^ у
k,"
у MpJ'
(x ) = fk,
\ jup I J у,/
■k," jup
dx dR[
dx dR[
dy
dRj
dy
где iup=i - 1 при Uj > 0, iup=i + 1 при u,,- < 0, jUp= j - 1 при v-j > 0, jup= j + 1 при v,j < 0. Полагается, что C30 =
0 при a10 = 1 и C03 =0 при a01 = 1, а также p10 = 0 при a10 = 0 и b01 = 0 при a01 = 0. При a10 = a01 = 0 функция Rk (x) является полиномом третьего порядка, в противном случае - рациональной функцией.
Аналогичные (6), (7) выражения применяются и для искомых функций, определенных на сетках S1 и S2. При этом значения компонент вектора скорости в узлах, где эти компоненты не определены, рассчитываются либо с использованием соответствующих интерполяционных функций вида (7), либо простым осреднением
v+1/2, j = (v" j-1/2 + V+1, j-1/2 + v"j+1/2 + v+1,-+1/2)/4, ui, j+1/2 = (ui-1/2, j + u,+1/2, j + u,-1/2, j+1 + ui+1/2,j+1)/4'
Как видно из (7), при введении сеточных
1 ч -/"к," -/"к," гк
функций jx j , jy j в дополнение к /, интерполяционную функцию Rjk(x) удается построить на шаблоне, состоящем лишь из четырех узлов. В противном случае для получения интерполяционной функции, по свойствам подобной Rj-^x), пришлось бы использовать шаблон из значительно большего числа точек, что менее удобно, особенно вблизи границ расчетной области. Следует отметить также, что (7) можно без изменений применять и на неравномерных сетках.
Еще одно преимущество дополнения набора искомых функций р, u, v, p, F системы (2) их пространственными производными , /к'" можно
продемонстрировать на следующем примере [14]. Пусть в одномерной задаче имеется волна с дли-
ч ч ч гк
ной, равной двум ячейкам сетки, такая, что = t = /,+1 = /,+2 , а знаки производных чередуются: s,g" (Лк;."1) = - s,g" () = s,g" (^ih) = -sig" (/к]"+1) . При использовании интерполяционной
функции, подобной по свойствам функции (7), но построенной с применением только узловых значе-
v ¡•к," гк¿~к,гг гк^ ^
ний -1 , /i , /,+1 , /,+2 самой интерполируемой функции, распознать наличие такой волны на отрезке x, < x < xi+1 не удастся, поскольку производные на концах этого отрезка окажутся равными
нулю. Чтобы распознать подобную волну при использовании такого подхода, понадобится более мелкая сетка. В то же время, если применяется интерполяционная функция (7), построенная по узловым значениям не только интерполируемой функции, но и ее производной, то подобная волна в решении на отрезке x, < x < xi+1 будет обнаружена. Это свойство интерполяционной функции (7) позволяет точнее описывать более резкие, в том числе и близкие к разрывным, изменения искомой функции на более грубых сетках.
Полулагранжевы методы решения уравнения переноса основаны на том, что значение функции f, удовлетворяющей уравнению
f + u-Vf = 0, остается вдоль траектории xa(t) произвольной частицы среды a постоянным: f (xa (t)) = const. Отсюда f (Xan) = f (Xa(tn)) = f (Xa(tn - Atn)). При достаточно малых Atn имеем xa(tn - Atn) » xa (tn) - Atn ua (tn) = xan
л-аК^ ¿-и j -л-a '
, так что можно записать
/(хЛ »/(X/ - Ы" иа"). На основе этого выражения и получены равенства (6). Учитывая, что применение равенств (6) связано с приближенной заменой ха(Г - А/п) на ха" - Ы" иа", а также имея в виду, что при определении интерполяционной функции Я/(х) согласно (7) значения производных /к;, /у,;" в узле (ц,, /р) не используются, шаг по времени Агп естественно ограничить условием
At" <-
1
í
max|u"|
— +—
Ax Ay
V
(10)
(11)
Одномерный аналог этого условия имеет вид A n Ax
Atn <-
max | u |
i
1.4. Решение уравнений неконвективной части
Уравнения ff = Gk неконвективной части (5) аппроксимируются по времени следующими неявными соотношениями
Р^1 -Р
At" u"+1 - u*
At"
p"+1 - P At"
= -p*V •u"+1, = -(p* )-1 Vp"+1, - = -p*C*2 V • u"+1.
Я+1 —*
(12)
Наряду с этим полагается ^ = ^ , поскольку неконвективная часть в уравнении для ^ отсутствует. При такой аппроксимации схема расчета р,;п+1, и;п+1, р"1 не имеет ограничений на шаг по времени. Входящие в (12) пространственные производные приближаются центральными разностями на разнесенных сетках 80-82.
Применив ко второму уравнению системы (12) операцию дивергенции и выразив У ип +1 из третьего уравнения, легко получить
„и+1
V^
( Vp"+1 ^
Р
p*c;2 At
."2
+ -
V • u*
At"
(13)
Пространственная дискретизация равенства (13) определяется пространственной дискретизацией уравнений (12).
Нетрудно заметить, что в случае несжимаемой жидкости, когда скорость звука С$ равна бесконечности, уравнение (13) сводится к эллиптическому уравнению У - (Ур"+1 / р*) = У - и* / АГ.
Расчет Pj
ijj , Pjj начинается с опреде-
ления р; из уравнения (13). При этом система линейных алгебраических уравнений относительно узловых значений давления р/^1, которая получается в результате пространственной дискретизации (13), решается итерационным методом последовательной верхней релаксации. Затем из второго уравнения (12) находится и/+1, а р;п+1 определяется из соотношения
п+1_ *
р;+1 = р; +р - р
*2 CS ,ij
Изложенная процедура решения неконвективной части (12) исходных уравнений (1) на основе уравнения (13) известна как Unified Procedure [17], а в сочетании с одним из методов группы CIP для расчета конвективной части - как С1Р-СиР (CIP-Combined Unified Procedure) [14]. В методах ICE и PISO, которые являются обобщением на сжимаемые среды методов расчета несжимаемых сред, для нахождения давления используются аналогичные (13) уравнения. Однако в отличие от них в методе CIP-CUP исходные уравнения записываются относительно неконсервативных переменных p, u, p и применяется их расщепление на конвективную и неконвективную части. В результате уравнение для давления (13), используемое в методе СIP-СUP, оказывается проще. Метод dP-GUP позволяет рассчитывать не только течения с широким диапазоном проявления сжимаемости среды, но и многофазные течения с большой разницей акустических импедансов фаз. Последнее обусловлено тем, что (13) обеспечивает непрерывность отношения Vp/p (которое, в соответствии со вторым уравнением из (12), можно рассматривать как ускорение) на границе раздела сред даже при отношении их плотностей порядка 103, в отличие от методов ISE и PISO, которые этого не гарантируют [14].
Ур авнения = Gkx - ux •Vfk,
Vfk неконвективной части (5) для
/к = о - и у
искомых производных, определяемых на основной сетке 80, аппроксимируются следующими соотношениями
- At" u "
a
f kn+1 f к * J x,j_J x,y
Atn
un - un vn - vn
_ г к * " i+1, j "i-1, j - г к * W+1, j Vi-1, j j x,ij ^ A J y,
кк Gi+1,j - Gi-1,j
2Ax
2Ax
' y,i
Gk.+, - G*. ,
i, j+1 i, j-1
2Ax
к n +1 к *
J у M_J yM
Atn
un - un vn - vn
г к * i, j+1 i, j-1 f к * i,j+1 4j-1 x,j ТА-.. У,И
(14)
2Ay
где
^к
2Ay y,j 2Ay
Gi±1, j =( f±j - fix*, j )/Atn , G?.±1 =(f-ji1 - fja)/Atn . Такая аппроксимация для Gi: основывается на том, что уравнение fk = Gk неконвективной части (5) аппроксимируется выражением (/к"+1 - fk*)/At" = Gk . Соотношения (14) применяются в конце алгоритма расчета
неконвективной стадии, когда рассчитаны все значения fjkn+1. Аналогичные (14) выражения применяются и для искомых производных, определяемых на сетках S1, S2. При этом как отмечалось выше, значения компонент вектора скорости в узлах, где эти компоненты не определены, рассчитываются либо с использованием соответствующих интерполяционных функций вида (7), либо простым осреднением вида (9).
1.5. Искусственная вязкость Схема CIP-CUP не является консервативной, поэтому при решении задач с интенсивными ударными волнами требуется введение искусственной вязкости, которое сводится к тому, что в правую часть второго и третьего уравнений системы (1) добавляются слагаемые Qu и Qp, соответственно
Vp
ut + u Vu = —— + Qu, P
p, + u-Vp = -pC2 Vu + Qp,
где
Qu =-1 Vgv, Qp =-(k-1)qv Vu^,
P p
f K +1
qv = с P I^-Cs AU +— AU2
AU = min (0,1V-u), (15)
к = фГ + (1-ф) g, Cv = ф cvJ + (1-ф) c,g ,
cv1, cv,g - коэффициенты искусственной вязкости жидкости и газа, l - параметр, имеющий порядок характерного размера ячеек расчетной сетки.
В результате такой коррекции исходных уравнений к описанным выше конвективной и неконвективной стадиям расчета добавляется стадия учета искусственной вязкости. Если полученные на неконвективной "невязкой" стадии значения скорости и давления обозначить u"+,L и p^+L (сеточные индексы i, j для краткости опущены), на стадии
учета искусственной вязкости будут применяться явные соотношения
u
= um + AtnQ*.
р-+1 = р:+± +А?ие;.
Сомножители О*, Q*p определяются через значения сеточных функций, рассчитанных на конвективной стадии. Коэффициенты искусственной вязкости обычно принимаются для газа в интервале (0.5, 1), а для жидкости - в интервале (0, 0.5).
1.6. Выбор шага по времени Уравнения (12) решаются неявным методом, не имеющим ограничения на шаг по времени, а ограничение (10) конвективной стадии является весьма слабым. Расчеты показывают, что и аппроксимация (14) уравнений /к = ОX -их -У/к , ¡к = Оку - и -У/к неконвективной части не накладывает каких-либо ограничений на устойчивость вычислений. Поэтому шаг по времени можно выбирать, исходя из особенностей рассматриваемых задач.
В задачах с крутыми (в том числе и ударными) волнами при использовании явных схем, например, схемы С. К. Годунова, шаг по времени обычно выбирают, опираясь на известное условие устойчивости Куранта. В частности, в случае равномерной эйлеровой сетки полагают
Atn =к
1
(
CRT max(| u". | +C"SJj)
11 — +—
Ax Ay
V
(16)
где КСЕТ - число Куранта, кСЕТ < 1. На практике обычно принимают КСЯТ ~ 1. Методика настоящей работы позволяет проводить расчеты и при КСДТ > 1 при условии выполнения (10). Однако при КСДТ > 1 может возникать чрезмерное размывание фронтов крутых волн. Одномерный аналог равенства (16) имеет вид
Ат
А(п = КС„Т-. (17)
СК1 /I п I , \ 4 у
тах(| и,. | +СБ1)
I '
Рис. 1 иллюстрирует зависимость погрешности §, возникающей при использовании методики настоящей работы для расчета ударных волн, от числа КСДТ в диапазоне 0.01 < КсДТ <10 в случае определения А1п согласно (17). Приведенные результаты получены при решении одномерной задачи о поршне, который в начальный момент мгновенно меняет свою скорость от 0 до 50 м/с и далее вдвигается с этой скоростью в покоящуюся жидкость (воду) ортогонально своей плоской поверхности. В жидкости возникает уходящая от поверхности поршня ударная волна. Параметры покоящейся жидкости: р =1 бар, р = 1000 кг/м3, С8 = 1482 м/с, Г = 7.15, В = 3.072 103 бар. Длина расчетной области составляет 3 м. В виду наличия в жидкости ударной волны используется искусственная
вязкость с коэффициентом су1 = 0.3. Погрешность § определяется выражением
§=№
I,Ягет г г,Ыит
)2
(18)
где р1Яг,ет - аналитическое решение (давление) и
р11иш - численное решение в г-м узле сетки, N -
число узлов. Приведенные на рисунке значения § соответствуют моменту времени / = 0.42 мс. Расчеты выполнены на сетках с числом узлов N = 200, 2000. В обоих вариантах условие (10) конвективной стадии не нарушается.
10"2 0.1 1 кгк, Ю Рис. 1. Зависимости погрешности § (18) методики настоящей работы при расчете задачи о поршне, вдвигаемом в покоящуюся жидкость с образованием в ней ударной волны, от числа КСЯТ (17). Сетки равномерные из 200 (кривая 1) и 2000 (кривая 2) ячеек.
На рис. 1 видно, что с увеличением КсЯТ погрешность счета § во всем диапазоне 0.01 < КСЯТ < 10 монотонно увеличивается. При этом оптимальное значение КСЯТ находится в интервале 0.1 < КСЯТ < 1, поскольку при КсЯТ < 0.1 погрешность практически не уменьшается, а при КСЯТ > 1 она с ростом КСЯТ довольно быстро увеличивается. Эти особенности учитывались при численном решении представленных ниже задач. 1.7. Дробление шага по времени на неконвективной стадии
Используемое в настоящей работе расщепление на конвективную и неконвективную стадии позволяет естественным образом проводить расчеты этих стадий с разными шагами по времени [28]. При решении некоторых задач это оказывается полезным, приводя к повышению экономичности вычислений, в частности, в тех случаях, когда на неконвективной стадии требуется существенно меньший шаг по времени, чем на конвективной стадии.
В операторном виде расчет одного временного шага по используемой в настоящей работе методике с единым для обеих стадий шагом А/п можно схематично представить как
Я* = Ц(А/п )яп, < = 12(АТ )я* , где я = (П fx, fy), f = (р, и, у, р, ^)Т и /.(А/*) , Ь2(А/п) - дискретные операторы соответственно для конвективной и неконвективной частей.
При использовании на неконвективной стадии дробного шага будем иметь
Я* = Ц (А/п)яп, яп+1 = [ц (А/п / К)]" я*, где К - положительное целое число. В этом случае после расчета конвективной части с шагом А/п неконвективная часть рассчитывается К раз с шагом А/п/К
В работе [28] шаг по времени А/п предлагается определять как
(
А/п = тт
И
И
"1 тах(| <|)
2 тах^.)
; )
(19)
где И - шаг равномерной сетки, с1, с2 - параметры с типичными значениями с1 = 0.2, с2 = 10 (в случае несущественного влияния сжимаемости или очень длинных волн с2 = да). Величина дробного шага неконвективной стадии определяется пространственной неоднородностью скорости звука вне контактной границы
( И
К
= тт
тах(| С" , - С" |)'
,1+1,; И
А/п
тах(| С
п ¡-т
в ,1, ;+1 -
(20)
)
При с3 = да конвективная и неконвективная стадии рассчитываются с одинаковым шагом.
Эмпирический подбор значений трех параметров с1 - с3 при использовании (19), (20) представляется более сложным, чем подбор одного параметра в случае применения (16) или (17). Поэтому в настоящей работе шаг по времени определяется, в основном, согласно (17) в одномерном случае и (16) в двумерном.
2. Тестовые расчеты
Работа алгоритма и программы расчета контролировалась путем сопоставления результатов численного решения ряда тестовых задач с аналитическими решениями или численными решениями других авторов. Результаты расчетов некоторых тестов приведены ниже.
2.1 Распространение возмущения квадратной формы
Для тестирования алгоритма решения уравнений переноса рассмотрена плоская задача о переносе профиля, изображенного на рис. 2а, с постоянной скоростью в диагональном направлении относительно координатных осей. Изначально профиль представляет собой возмущение с единичной амплитудой нулевого значения функции ф в области квадратной формы. Распространение возмущения описывается уравнением переноса // + и У^ = 0 , где Р = ф, либо ^ = гя[(1 -е)р(ф-0.5)], е = 0.01. Задача решается в безразмерных переменных. При /
с
3
с
= 0: ф = 1 в квадратной области 20x20 ячеек в центре расчетной области размером 160x160 ячеек, ф = 0 в остальной ее части (рис. 2а). Сетка равномерная с Ах = Ду = 1, шаг по времени Дг = 0.2, скорость и = (1, 1).
В рамках контактного взаимодействия эту задачу можно интерпретировать как движение в поле постоянной скорости несжимаемой среды с квадратным включением, занятым другой несжимаемой средой, при отсутствии вязкости и поверхностного натяжения. Значения 0 и 1 функции ф идентифицируют области, занятые разными средами. Граница между средами находится там, где 0 < ф < 1.
Рис. 2. Распространение квадратного возмущения. Профиль и изолинии функции ф при t = 0 (а) и t = 300Дí (б-
е). Результаты получены при Р = ф (б, в); Р = гя [(1 -е) Я ( ф-0.5) ] (г-е) с использованием интерполяционной функции (7) при а10 = а01 = 0 (полином третьего порядка) (б); интерполяционной функции (7) при а10, а01 согласно (8) (в, г ); схемы Лакса-Вендроффа (д); схемы «против потока» (е).
Согласно начальной дискретизации функция ф изменяется от 0 до 1 в области шириной в одну ячейку. Поэтому в поле изолиний ф начальная граница возмущенной области (начальная граница между средами) выглядит как четкий квадратный контур (рис. 2а). На рис. 2б приведены результаты, полученные без тангенциального преобразования (Р = ф) с использованием интерполяционной функции (7) при а10 = а01 = 0, являющейся в этом случае полиномом третьего порядка. Видно, что профиль ф несколько размывается (так что неопределенность положения границы между средами увеличивается), но его форма остается близкой к квадратной. Также возникают небольшие осцилляции в окрестности границы возмущения, где градиент ф испытывает наибольшие изменения. При использовании интерполяционной функции (7) с определением а10, а01 согласно (8), то есть с возможностью переключения между кубическим полиномом и рациональной функцией, решение становится монотонным (рис. 2в). Это объясняется тем, что при использовании рациональной функции - отношения полиномов - порядок аппроксимации уравнения переноса понижается. При дополнительном применении тангенциального преобразования
Р = гя [(1 -е) я (ф-0.5)] профиль функции ф через
300 временных шагов почти полностью совпадает с начальным, и граница возмущенной области (граница между средами) остается четкой (рис. 2г). Нужно отметить, что не для всех схем тангенциальное преобразование приводит к повышению эффективности расчета уравнения переноса. Это иллюстрируют рис. 2д, е с результатами расчета по схемам Лакса-Вендроффа и «против потока» в сочетании с тангенциальным преобразованием.
Таким образом, используемый в настоящей работе метод расчета конвективных слагаемых характеризуется как относительно небольшой численной диффузией, так и монотонностью решения в областях его больших градиентов и их быстрого изменения. Эти качества - результат сочетания по-лулагранжева подхода с интерполяционной функцией (7), принимающей, в зависимости от локальных особенностей решения, форму либо кубического полинома, либо рациональной функции. При дополнительном использовании тангенциального преобразования функции-идентификатора ф в задачах контактного взаимодействия граница между контактирующими средами, определяемая ненулевым градиентом ф, размывается не более чем на одну-две ячейки.
2.2 Взаимодействие акустической волны с контактной границей
Для тестирования алгоритма расчета уравнений неконвективной части (5) решалась одномерная задача о взаимодействии акустической волны с границей раздела сред типа газ-жидкость [28]. Задача решалась в безразмерных переменных. Область газа: 0 < х < 0.5, область жидкости 0.5 < х < 1. Для газа у = 1.4, для жидкости Г = 7.15, В = 3072.04. В начальный момент времени газ и жидкость покоятся (и = 0) при давлении р =1, плотность газа р = 1.025*10-3, плотность жидкости р = 1. Акустическая волна входит в область газа через ее левую границу х = 0, где задается р = 1 + 0.1(1 - 008 Юt),
Ю= 2яС°,8 /1, С° » 37 - невозмущенная скорость звука в газе, 1 = 0.225 - длина волны. Расчеты проводились при Дх = 1.25 • 10-3, Дt = 5 • 10-6 на сетке из 800 ячеек. Как и в работе [28], вычисления выполнены без учета конвективных слагаемых.
На рис. 3 видно, что акустическая волна при распространении по газу практически не размывается. В процессе взаимодействия с границей раздела между газом и жидкостью она частично отражается, а частично проходит через нее. Результаты настоящей работы хорошо согласуются с результатами работы [28].
2.3. Распад разрыва
Для тестирования алгоритма расчета в целом, т.е. совместной работы его конвективной и неконвективной составляющих, решались три одномер-
ные задачи о распаде разрыва. В первой задаче в начальном распределении параметров однородной среды (жидкости) имеется разрыв скорости. Во второй и третьей задачах на границе контакта газ-жидкость терпят разрыв скорость и давление соответственно.
ширине размывания разрывов к-сят = 0.4 и су1 =0 (рис. 4в).
близки к случаю
О 0.5 х 1
Рис. 3. Взаимодействие акустической волны с границей контакта газа и жидкости (пунктирная линия). Распределения давления до взаимодействия в момент / »1.2 • 10-2
(штриховая линия), и в ходе взаимодействия в момент ? » 1.6 •Ю-2 (сплошная линия). Символы • - результаты работы [28].
2.3.1 Распад разрыва в однородной жидкости
В начальный момент в жидкости с плотностью р = 103 кг/м3 и давлением р =1 бар имеется разрыв скорости: и = и0 < 0 при х > х0 и и = 0 при х < х0. Эту одномерную задачу о распаде разрыва в однородной жидкости можно интерпретировать как задачу об ударе столба жидкости бесконечного радиуса по поверхности покоящейся жидкости с идентичными свойствами. Параметры жидкости: Г = 7.15, В = 3072 бар. Полагается х0 = 20 м, и0 = -50 и -500 м/с. Расчетная область длиной 40 м (по 20 м для подобластей с подвижной и покоящейся жидкостью) покрывается равномерной сеткой из 300 ячеек. Расчеты проводились при ксят = 0.1, 0.4 и 0.8 в (17) и су1 = 0, 0.2.
Результаты расчетов данной задачи представлены на рис. 4. Видно, что при относительно слабом ударе, когда и0 = -50 м/с (рис. 4а-г), применение искусственной вязкости не является необходимым. В этом случае при ксят = 0.1 и су,1 = 0 (рис. 4а) на фронте ударных волн, распространяющихся в обе стороны от разрыва, возникают небольшие осцилляции. Давление за ударными волнами и положение их фронтов хорошо совпадают с аналитическим решением. При увеличении ксят до 0.4 ширина размывания фронтов ударных волн увеличивается, а осцилляции на них отсутствуют (рис. 46). Дальнейшее увеличение ксят до 0.8 приводит к еще большему размыванию фронтов ударных волн (рис. 4в). Таким образом, в используемом в настоящей работе методе схемная вязкость с ростом ксят монотонно увеличивается. Введение искусственной вязкости также приводит к увеличению ширины размывания разрывов и уменьшению осцил-ляций на них. На рис. 4г представлены результаты, полученные при ксят = 0.1 и су1 = 0.2, которые по
400-, о,
се Ю
0' 400
а
Кскг= 0 1
СК1 0
1 1 1
5000-
0400-
[ \ 6
\ Кскт= 0.4 Су,/ = 0 \
■ 1 1
! \ ; 1 в
/ Кда^О.8 С>;/ 0 \
е
- . Ксяг= 0.4 ,
0
—■-1-^ 1
5000-
_ ж
• к0.8 *
£у.г о :
1 , , -г
Рис. 4. Пространственные распределения давления в момент í = 8 мс в задаче о распаде разрыва в однородной жидкости при и0 = -50 м/с (а-г), -500 м/с (д-з). Сплошные жирные линии - аналитическое решение, тонкие линии с символами • и штриховая линия - численные решения настоящей работы.
При увеличении скорости движущейся жидкости до и0 = - 500 м/с (рис. 4д-з) образуются более интенсивные ударные волны. В таком случае одной схемной вязкости, которую можно увеличить посредством увеличения ксят (т.е. путем увеличения шага по времени), уже недостаточно для получения удовлетворительного численного решения. Как видно на рис. 4д-ж, без введения искусственной вязкости давление за ударными волнами и положение их фронтов рассчитываются с большой погрешностью. Увеличение числа ксят при сУ1 = 0 способствует большему размыванию фронтов ударных волн и уменьшению осцилляций в их окрестности. В результате введения искусственной вязкости с су1 = 0.2 в случае ксят = 0.1 давление за фронтом ударных волн и положение самих фронтов становятся довольно близкими к аналитическим (рис. 4з). Для ксят = 0.8 искусственной вязкости при сУ1 =0.2 недостаточно (рис. 4з): численное решение ближе к аналитическому, чем при су,1 =0 (рис. 4ж), но все еще заметно отличается от него. Для ксят =0.8 удовлетворительное согласование давления за ударными волнами и положения их фронтов с аналитическим решением можно получить при су,1 = 1.
2.3.2 Распад разрыва, обусловленного скачком скорости на границе газ-жидкость
Рассматривается одномерная задача о распаде разрыва на границе газ-жидкость, в которой при t = 0 давление в обеих средах одинаково (р = 1 бар), а скорость терпит разрыв. В начальный момент времени положение контактной границы х0 =35 м, в области жидкости (х > х0) и = и0 < 0, р = 103 кг/м3, ф = 1, Г = 7.15, В = 3072 бар, в области газа (х < х0) и = 0, р =1 кг/м3, ф = 0, у = 1.4. Данную задачу можно интерпретировать также как задачу об ударе столба жидкости бесконечного радиуса по покоящемуся газу. Расчеты проводятся с использованием преобразования Р = гя[(1 -е)Я(ф- 0.5)], с искусственной вязкостью при Су,1 = 0.2, = 0.6. Расчетная область длиной 70 м (по 35 м на подобласти с газом и жидкостью) покрывалась равномерной сеткой из 300 ячеек.
Результаты расчетов представлены на рис. 5. Взаимодействие движущейся жидкости с неподвижным газом приводит к тому, что в газе и жидкости образуются ударные волны, контактная граница смещается в сторону газа. При этом в расчетах с к-сят = 0.2 как при и0 = -50 м/с (рис. 5а), так и при и0 = -500 м/с ( 5г) возникают большие отклонения давления от аналитического решения в окрестности фронта ударной волны в жидкости. Так, в момент ^ = 4 мс при и0 = -50 м/с величина отклонения примерно в 3 раза превышает приращение давления на фронте ударной волны в аналитическом решении, а при и0 = -500 м/с - уже в 7 раз. Как видно, по мере распространения ударной волны в жидкости отклонение от аналитического решения уменьшается (момент Возникновения подобных отклонений не удается избежать при варьировании вязкости, как схемной (изменение ксят), так и искусственной (изменение сУ1, су>г). Интересно отметить, что при расчете без учета конвекции параметров, характеризующих свойства среды (плотности и функции-идентификатора), таких больших отклонений не возникает. Поэтому можно предположить, что в процессе формирования в численном решении ударных волн в малой окрестности контактной границы размывание профиля плотности на контакте, возникающее при расчете конвекции, приводит к появлению сильных возмущений давления.
Чтобы избежать этого, можно попытаться воспользоваться описанной в пункте 1.7 техникой дробления временного шага. При этом расчет одного шага ДГ конвективной части предваряется расчетом К шагов неконвективной части с дробным шагом ДЛК. Результаты применения подобной техники приведены для и0 = -50 м/с на рис. 5б и для и0 = -500 м/с на рис. 5д. В обоих случаях на первом шаге дробление не применяется. Единый для конвективной и неконвективной частей первый шаг Д^ определяется согласно (17) при ксят = 0.2. В случае и0 = -50 м/с, начиная со второго шага и
далее, используются (19), (20) с коэффициентами с1 = 0.4, с2 = 10, с3 = 0.1. Получающийся суммарный шаг Л/2 примерно в 50 раз превышает шаг, определяемый условием (17) при Kcrt = 0.2. Сначала неконвективная часть рассчитывается около 80 раз с дробным шагом ~Л/2/80. В результате в численном решении формируются расходящиеся от контактной границы ударные волны, и эти волны успевают отойти от нее на расстояние, примерно равное 2 ячейкам в газовой среде и 10 ячейкам в жидкой среде. Затем рассчитывается конвективная часть с суммарным шагом Л/2, который в данном случае удовлетворяет условию (11). Третий, и последующие суммарные шаги от шага Л/2 существенно не отличаются, а количество их разбиений для расчета неконвективной части уменьшается. В приведенные на рис. 5б моменты времени /1 и /2 это количество соответственно равно 50 и 10. Как видно на рис. 5б, при использовании такого подхода полностью избежать отклонений от аналитического решения не удается. Однако их величина существенно уменьшается по сравнению с тем, что представлено на рис. 5а. В случае u0 = -500 м/с заметно снизить отклонения от аналитического решения на втором расчетном шаге удается лишь при с1 » 10, с2 » 10, с3 » 0.08. Суммарный шаг расчета конвективной части Л/2 примерно в 70 раз превышает величину шага, определяемого условием (17) при kCRT = 0.2. При расчете неконвективной части он делится на 100 дробных шагов. Как и в случае u0 = -50 м/с, при расчете только неконвективной части за время Л/2 формирующиеся ударные волны в газе и жидкости отходят от контактной границы на расстояние порядка 2 и 10 ячеек, соответственно. Однако при такой величине шага Л/2 ограничение (11) на величину временного шага конвективной стадии нарушается. Многократное нарушение условия (11) может привести к большой погрешности, поэтому, начиная с третьего шага, расчет варианта u0 = -500 м/с продолжается по обычной схеме с единым для конвективной и неконвективной стадии шагом по времени, определяемым согласно (17) при kCRT = 0.2. Как видно на рис. 5д, несмотря на имеющееся однократное нарушение условия (11), отклонения давления от аналитического решения для u0 = -500 м/с оказываются по сравнению с рис. 5г существенно меньшими.
Для уменьшения больших погрешностей численного решения в данной задаче можно также воспользоваться коррекцией плотности, целью которой является уменьшение размывания профиля плотности на контактной границе при расчете конвекции. Суть коррекции заключается в том, что после расчета конвекции значения плотности p, и ее производной px¡ в узле, прилегающем к контактной границе с «подветренной» стороны (в котором выполняется условие (ф,- - 0.5)-(фщР - 0.5) <
0), заменяются на значения р*корр и рхк°рр, определяемые следующим образом р;корр = рП, р;т = 0 при (ф* - 0.5) -(ф -0.5) > 0 ,
р*= рПр, рхг = 0 при (ф; -0.5) •(фП -0.5) < 0,
где индексами п и * обозначены величины соответственно до и после расчета конвективной части. Такой способ предотвращения размывания профиля плотности на контактной границе более чем на одну ячейку является довольно простым (по сравнению, например, со способом, предложенным в [34]). Несмотря на это, в рассматриваемой задаче он позволяет уменьшить отклонения от аналитического решения так, что при обоих значениях и0 в момент ^ они становятся уже практически незаметными (рис. 5 в, е).
Рис. 5. Пространственные распределения давления при t = t2 и плотности при t = t2 в задаче о распаде разрыва на границе газ-жидкость, ^ = 4 мс, t2 = 20 мс. При t = 0 в области газа (х < х0) и = 0, в области жидкости (х > х0) и = -50 м/с (а-в) и -500 м/с (г-е). Рисунки а, г - шаг по времени определяется по формуле (17) при ксят = 0.2; б, д -с дроблением шага согласно (19), (20) при с1 = 0.4, с2 = 10, с3 = 0.1 (б) и с1 = 10, с2 = 10, с3 = 0.08 (д); в, е - с коррекцией плотности, и использованием (17) при ксят = 0.2.
Сплошные кривые - аналитическое решение, тонкие кривые с символами • - численное решение настоящей работы, 1 - контактная граница, 2 - ударная волна.
Следует отметить, что в задаче распада разрыва, обусловленного скачком скорости на границе газ-жидкость оба описанных приема сокращения погрешности решения требуют подбора оптимального сочетания значений некоторых параметров. При использовании дробления временного шага на неконвективной части - это параметры с1-с3 и коэффициенты искусственной вязкости су1, са при использовании коррекции плотности - это ксят, с„/, с„8. 2.3.3 Распад разрыва, обусловленного скачком давления на границе газ-жидкость Рассматривается одномерная задача о распаде разрыва на границе газ-жидкость, в которой при t =
0 обе среды покоятся (и = 0), а давление терпит разрыв. В начальный момент времени положение контактной границы х0 = 360 м, в области жидкости (х > х0) р = 14087.9 бар, р = 1216 кг/м3, ф = 1, Г = 7.15, В = 3072 бар, в области газа (х < х0) р = 1 бар, р = 1 кг/м3, ф = 0, у = 1.4. Расчетная область длиной 1700 м (360 м - подобласть с газом и 1340 м - подобласть с жидкостью) покрывалась равномерной сеткой из 300 ячеек. Расчеты проводились при ксят = 0.5 и су1 = 0, с^ = 0.7 (без использования коррекции плотности и дробления шага на неконвективной части). Результаты расчетов представлены на рис. 6.
Рис. 6. Пространственные распределения давления и плотности при t = 0.4 с в задаче о распаде разрыва на границе газ-жидкость. При t = 0 в области газа (х < х0) р =
1 бар, в области жидкости (х > х0) р = 14087.9 бар. Сплошные кривые - аналитическое решение, символы • - численное решение при = 0.7, штриховая кривая -численное решение при = 0 (без искусственной вязкости). Цифрами отмечены: 1 - контактная граница; 2 -волна разрежения в жидкости; 3 - ударная волна в газе.
Как видно, при использовании искусственной вязкости в газе согласование численного решения настоящей работы с аналитическим оказывается вполне удовлетворительным и в области газа, и в области жидкости, и в окрестности контактной границы. Если же искусственную вязкость не использовать, то давление в окрестности фронта ударной волны в газе заметно отличается от аналитического. Положение фронта ударной волны в этом случае также определяется с существенной погрешностью.
2.4. Цилиндрический взрыв
Для тестирования двумерного алгоритма рассматривались двумерные задачи о цилиндрическом взрыве с учетом и без учета возникающей при взрыве неоднородности среды [24]. 2.4.1 Цилиндрический взрыв без учета неоднородности среды
Рассматривается цилиндрический взрыв в покоящемся газе: и = 0, у = 1.4. Задача решается в безразмерных переменных. При t = 0: р = 1, р = 3. в
круговой области г (х - 0.5)2 + (у - 0.5)2 < 0.08 и р = 0.1, р = 0.1. в области г > 0.08. При t > 0 в окрестности границы г = 0.08 области взрыва формируются радиально сходящаяся к точке г = 0 (х = 0.5, у = 0.5) волна разрежения и радиально расходящаяся от этой точки ударная волна. Расчетная область (0
< х < 1)х(0 < у < 1) покрывается сеткой из 100x100 ячеек с Ах = Ау = 0.01. Шаг по времени = 5 • 10-4. Результаты расчетов представлены на рис. 7. Начальная граница области взрыва (пунктирная кривая 1 на рис. 7а) представляет собой аппроксимацию окружности на равномерной декартовой сетке отрезками прямых линий.
На рис. 76 видно, что при использовании искусственной вязкости (15) с коэффициентом су^ = 1.4 численное решение в сечении у = 0.5 получается довольно близким к результатам работы [24], которые были получены с применением явной схемы для расчета неконвективной части. Изолинии плотности в области фронта ударной волны и в окрестности заднего фронта волны разрежения в численном решении настоящей работы получаются близкими к круговым (рис. 7а). Относительно небольшое различие радиальных профилей плотности в сечениях у = 0.5 и х = у обусловлено тем, что расчетная сетка является недостаточно мелкой. При этом в обоих сечениях ширина размывания фронта ударной волны практически одинакова.
Для двумерных расчетов часто используют вариант искусственной вязкости, основанный на расщеплении по координатным направлениям
Qu = (Qu,Qv), Qu = -,Q p dx
du + dv
dx ^ dy
Q„ = -(K-1)| qvx
q.
= c,x p i -CSDu + (Du)2
j P dy '
% = c,y p| -CsAv + -K+1 (Av)2
Au = min(0, Idu/dx), Av = min(0, Idv/dy) (21)
На рис. 7в, г представлены результаты применения искусственной вязкости с пространственным расщеплением (21) при cvx = cvy = 1.4. Различие радиальных профилей плотности, в том числе положения фронта расходящейся ударной волны и ши-
рины его размывания, в сечениях у = 0.5 и х = у становится заметно больше, чем при искусственной вязкости без расщепления (15). При этом форма фронта ударной волны довольно сильно отличается от круговой. Указанные отличия свидетельствуют о том, что при использовании вязкости вида (21) ее влияние в диагональном направлении оказывается слабее, чем в направлении координатных осей.
Если искусственную вязкость не применять вообще (рис. 7д, е), то ширина размывания фронта расходящейся ударной волны становится значительно меньше. Однако при этом в окрестности фронта возникают довольно большие осцилляции.
2.4.2 Цилиндрический взрыв с учетом неоднородности среды
Постановка задачи и входные данные описаны в предыдущем параграфе. Отличие состоит в том, что при / = 0 в круговой области взрыва принимается у = 1.8, ф = 0, а вне этой области у = 1.2, ф = 1. Чтобы приблизить условия расчета к работе [24], тангенциальное преобразование функции ф не использовалось, то есть расчеты проводились при ^ = ф. Результаты представлены на рис. 8.
Рис. 8а-в получены при использовании двумерного варианта искусственной вязкости (15). На рис. 8в видно, что рассчитанный в настоящей работе радиальный профиль плотности в сечении у = 0.5 хорошо согласуется с профилем, полученным в работе [24]. Согласно рис. 8а ширина размывания фронта ударной волны примерно одинакова в диагональном сечении х = у и в сечении у = 0.5. Контактная граница перемещается в обоих сечениях с близкой скоростью, так что ее форма сохраняется подобной начальной границе, которая представляет собой аппроксимацию окружности отрезками прямых (слошная кривая 1).
вой-верхней четверти расчетной области; нижний ряд - радиальные распределения плотности. Рисунки а, 6 - с искусственной вязкостью (15); в, г - с расщеплением искусственной вязкости по координатным направлениям (21); д, е - без искусственной вязкости. Пунктирная кривая 1 - начальная граница области взрыва (г = 0.08). Радиальные профили плотности в нижнем ряду соответствуют двум сечениям фрагментов расчетной области из верхнего ряда: у = 0.5 (символы •) и х = у (символы х), жирная сплошная кривая (6) - результаты работы [24] в сечении у = 0.5.
-•-1 ---1 -1-1
О Г 0.3 О Г 0.3 0 Г 0,3
Рис. 8. Задача о цилиндрическом взрыве с учетом неоднородности среды, момент t = 0.08. Верхний ряд - изолинии плотности в правой-верхней четверти расчетной области; средний ряд - начальное (кривая 1) и текущее положение контактной границы, определяемое изолиниями ф = 0.5; нижний ряд - радиальные распределения плотности. Рисунки а-в - с искусственной вязкостью (15); г-е - с расщеплением искусственной вязкости по координатным направлениям (21); ж-и - без искусственной вязкости. Обозначения те же, что и на рис. 7.
Применение искусственной вязкости в виде (21) с расщеплением по координатным направлениям (рис. 8г-е) в случае неоднородной среды приводит к изменениям в расчете ударной волны, во многом аналогичным случаю однородной среды (рис. 7в, г). Круговая форма контактной границы искажается, как и форма ударной волны. Это обусловлено тем, что контактная граница в диагональном направлении движется медленнее, чем вдоль координатных осей.
Без применения искусственной вязкости (рис. 8ж-и) фронт ударной волны сохраняется резким, но в области между контактной границей и ударной волной возникают осцилляции с довольно большой частотой и амплитудой. Форма контактной границы по сравнению с начальной формой становится более гладкой и близкой к окружности.
3. Иллюстрация применения методики
Для иллюстрации возможностей разработанной методики рассматривается задача об ударном воздействии высокоскоростной струи жидкости по твердой стенке. Предполагается, что параметры струи соответствуют параметрам кумулятивной струйки жидкости, возникающей на поверхности кавитационного пузырька при его схлопывании вблизи стенки. Изучение воздействия таких струек на стенку представляет значительный интерес в
связи с проблемой кавитационного разрушения поверхностей.
Расчеты [35] показывают, что наибольшие скорости струйки достигаются в случае, когда ее диаметр мал по сравнению с размерами пузырька (рис. 9а). Поэтому на начальном этапе исследования влиянием верхней и боковых границ пузырька можно пренебречь. В результате ударное воздействие струйки можно моделировать как удар столбика жидкости по стенке (рис. 9б).
Рис. 9. Схема ударного воздействия пузырька на стенку.
Рассматривается осесимметричная задача об ортогональном ударе цилиндрической струи жидкости с полусферическим концом по неподвижной плоской твердой стенке (рис. 9б). Струя считается бесконечно длинной и окруженной неограниченным объемом газа. Эффекты вязкости, теплопроводности жидкости и газа, а также влияние поверхностного натяжения не учитываются. Для жидкости началь-
ные значения параметров полагаются следующими: p = 103 кг/м3, Г = 7.15, B = 3047 бар, j = 1; а для газа - следующими: p = 1 кг/м3, g = 1.4, j = 0. В начальный момент времени давление всюду равно 1 бар. Радиус струи R = 25 мкм, ее скорость 500 м/с. Эти значения соответствуют кумулятивной струйке жидкости, образующейся на поверхности кавитационно-го пузырька в воде при комнатных условиях в том случае, когда пузырек примыкает к телу и в начале схлопывания имеет вид эллипсоида вращения с отношением полуосей ~0.84. Коэффициент искусственной вязкости в жидкости cvj = 0.4, а в газе cvg = 0.9, шаг по времени определяется из (16) при KCRT = 0.9. Расчетная область (0 < r/R < 0.5)х(0 < z/R < 0.25) покрывается немного сгущающейся к стенке расчетной сеткой 640x240 ячеек.
Одна из важных особенностей рассматриваемой задачи заключается в том, что натекание струи с выпуклым концом на стенку сменяется в некоторый момент ее растеканием с образованием тонкого пристеночного радиально разлетающегося с большой скоростью слоя жидкости (пристеночной струйки) [36, 37]. В результате этого граница струи подвергается сильным деформациям. Использование адаптивно-подвижных сеток, оптимальных в задачах с подвижными умеренно деформирующимися границами, в этой задаче становится проблематичным. Это связано с тем, что в конце стадии натекания струи расчетная сетка сильно искажается, требуется ее неоднократная перестройка, в том числе с изменением топологии, и последующая интерполяции решения со старой сетки на новую. В таких условиях эффективность применения адаптивно-подвижных сеток существенно снижается.
На рис. 10а, б показано сопоставление расчетных сеток, применяемых в настоящей работе (рис. 10а) и в расчетах работы [36] с применением адаптивно-подвижных сеток (рис. 10б), а также формы струи в конце стадии ее натекания, за которой следует растекание. В работе [36] наличие газа, окружающего струю, не учитывалось, то есть граница струи являлась границей расчетной области. На рис. 10в, г дано сопоставление изолиний давления и полей скорости. Как видно, положение границы струи (кривая 3), а также положение и криволинейная форма ударной волны в струе (кривая 1), полученные в настоящей работе, хорошо согласуются с результатами работы [36]. При этом важно отметить, что в настоящей работе начальный полукруглый профиль струи аппроксимируется на эйлеровой сетке отрезками прямых. Естественно, что вносимые таким образом погрешности зависят от размеров ячеек и по мере измельчения сетки уменьшаются.
В работе [36] была рассчитана только стадия натекания струи. Расчет стадии растекания не осуществлялся в связи со значительными сложностями перестройки сеток. В то же время методика настоящей работы позволяет рассчитывать не только
всю стадию натекания струи, но и ее последующее растекание с формированием пристеночной струйки (рис. 11). При этом каких-либо дополнительных трудностей не возникает.
Рис. 10. Разреженные в 8 раз расчетные сетки в малой окрестности области взаимодействия струи с полусферическим концом со стенкой в конце стадии натекания: а -
эйлерова сетка настоящей работы; 6 - адаптивно-подвижная сетка работы [36]. Жирная кривая - граница струи. Прореженное поле скорости, изолинии давления и на рисунке в изолинии функции-идентификатора: в -расчет настоящей работы; г - результаты работы [36]. Цифрами отмечены: 1 - ударная волна; 2 - волна разрежения; 3 - граница струи.
На рис. 12а приведено сравнение полученных в настоящей работе и в работе [36] распределений давления на стенке. Представлены четыре последовательных момента времени, включая момент растекания струи (которому соответствует рис. 11), рассчитанный только в настоящей работе. Можно отметить хорошее согласование с результатами работы [36] как величины давления за фронтом ударной волны (кривые /1), так и распределений давления на стенке в моменты /1-/3. Некоторая немонотонность кривых давления, полученных в настоящей работе, а также меньшее по сравнению с [36] максимальное значение давления на пике в окрестности контакта границы струи со стенкой в момент /3 могут быть обусловлены аппроксимацией этой границы отрезками прямых линий и недостаточным измельчением сетки.
Можно заключить, что, несмотря на погрешности, связанные с аппроксимацией гладких кривых на эйлеровой сетке отрезками прямых линий, согласование результатов расчетов настоящей работы с расчетами работы [36] на стадии натекания струи вполне удовлетворительное.
Рис. 11. Растекание струи с полусферическим концом: 1 -ударная волна; 2 - волна разрежения; 3 - граница струи.
Рис. 12. Радиальные распределения давления на стенке в последовательные моменты времени t1 - t4 в случае удара струи с полусферическим концом: символы • - расчет настоящей работы; сплошные кривые - результаты из работы [36]. Моменту tз в конце стадии натекания струи соответствует рис. 10, моменту t4 на стадии растекания соответствует рис. 11.
Заключение
Дан краткий обзор работ по численному моделированию задач контактного взаимодействия сжимаемых сред с большой разницей акустических импедансов, в том числе при сильных деформациях контактной границы и при наличии в средах ударных волн.
Приведены основные положения разработанной методики расчета подобных задач на основе метода С1Р-СиР, в котором применяются прямоугольные эйлеровы сетки без явного выделения контактных границ. Исходные уравнения записываются относительно неконсервативных гидродинамических переменных (плотности р, скорости и и давления р) и функции-идентификатора среды ф. Применение этой функции дает возможность неявного отслеживания контактных границ. Такой подход позволяет избежать проблем, возникающих при использовании криволинейных адаптивно-подвижных сеток, (подстраиваемых под геометрию контактных границ), в частности, проблем, связанных с сильными деформациями границ и изменением их топологии. В методе расчета уравнение переноса для ф заменяется на аналогичное уравнение относительно функции Р - тангенциального преобразования от ф. Важной особенностью метода С1Р-СИР является расщепление расчета одного временного шага на несколько этапов. В случае взаимодействия невязких нетеплопроводных сред таких этапов три. На первом этапе учитываются конвективные слагаемые, на втором - оставшиеся неконвективные слагаемые. Третий этап вводится для учета искусственной вязкости, которая нужна при расчете интенсивных ударных волн, поскольку решаемые уравнения не имеют формы сохранения. Отличительная черта метода С1Р-СИР - дополнение набора неизвестных р, и, р и Р их пространственными производными по каждому из координатных направлений. Эти производные также рассматриваются как неизвестные функции и определяются решением уравнений, полученных дифференцированием уравнений для р, и, р и Р по соответствующей координате. В рамках используемого в подходе СТР-СИР полулагранжева метода
расчета конвекции дополнение вектора неизвестных задачи пространственными производными от искомых функций обеспечивает возможность построения интерполяционной функции третьего порядка на компактном шаблоне, состоящем из двух узлов в одномерном случае и четырех узлов в двумерном. Кроме того, такое дополнение дает возможность описывать возмущения с длинами волн порядка двух ячеек сетки.
Для контроля правильности разработанных методики и программы расчета в целом и их основных составляющих решен ряд тестовых задач, имеющих либо аналитическое решение, либо численное решение, полученное другими авторами, которые принимались в качестве эталонных решений. Это задачи о распространении возмущения квадратной формы; о взаимодействии акустической волны с контактной границей; о распаде разрыва в однородной жидкости; о распаде разрыва на границе газ-жидкость, когда разрыв обусловлен либо скачком скорости, либо скачком давления; о цилиндрическом взрыве как без учета, так и с учетом неоднородности среды. Показано, что во всех случаях, используя имеющуюся в методике свободу выбора конкретного алгоритма расчета (например, изменяя коэффициент искусственной вязкости, применяя дробление временного шага на отдельных этапах алгоритма), можно достичь удовлетворительного согласования с эталонными решениями.
Возможности разработанной методики проиллюстрированы на задаче об ударе осесимметричной высокоскоростной струи жидкости с полусферическим концом по жесткой стенке. Полученные результаты сопоставлены с известными численными решениями, рассчитанными методом адаптивно-подвижных сеток с явным выделением межфазной границы. Сопоставление проведено лишь на начальной стадии удара, поскольку дальнейшее применение метода адаптивно-подвижных сеток становилось проблематичным, в отличие от методики настоящей работы. Показано удовлетворительное согласование результатов.
Работа выполнена в рамках программы РАН и при поддержке РФФИ (проект 12-01-00341-a).
ЛИТЕРАТУРА
1. Blake J. R., Gibson D. C. Cavitation bubbles near boundaries // Annu. Rev. Fluid Mech. 1987. V. 19. P. 99-123.
2. Zhu S., Cocks F. H., Preminger G. M., Zhong P. The role of stress waves and cavitation in stone comminution in shock wave lithotripsy // Ultrasound Med. Biol. 2002. V. 28. Р. 661-671.
3. Benjamin T. B., Ellis A. T. The collapse of cavitation bubbles and the pressures thereby produced against solid boundaries // Phil. Trans. R. Soc. Lond. 1966. V. 260. №1110. P. 221-240.
4. Crum L. A. Surface oscillations and jet development in pulsating bubbles // J. de Physique. 1979. V. 40. N. 8. P. C8-285-C8-288.
5. Ohl C.-D., Kurz T., Geisler R. , Lindau O., Lauterborn W. Bubble dynamics, shock waves and sonoluminescence // Phil. Trans. R. Soc. Lond. 1999. V. 357. №1751. P. 269-294.
6. Philipp A., Lauterborn W. Cavitation erosion by single laser-produced bubbles // J. Fluid Mechanics. 1998. V. 361. P. 75-116.
7. Nourgaliev R. R., Dinh T. N., Theofanous T .G. Direct numerical simulation of compressible multiphase flows: interaction
of shock waves with dispersed multimaterial media // CD-ROM Proc. of the 5th International Conference on Multiphase Flow, ICMF'04, 2004. Paper N.494.
8. Бураго Н. Г., Кукуджанов В. Н. Обзор контактных алгоритмов // Известия РАН. МТТ. 2005. №1. C. 45-87.
9. Kadioglu S. Y., Sussman M., Osher S., Wright J. P., Kang M. A second order primitive preconditioner for solving all speed multiphase flows // J. Comput. Phys. 2005. V. 209. №2. P. 477-503.
10. Kim J., Lowengrub J. Interfaces and multicomponent fluids // Encycliopedia of Mathematical physics. Elsevier, 2006. P. 135-144.
11. Donghong W., Ning Z., Ou H., Jianming L. A Ghost Fluid based Front Tracking Method for multimedium compressible flows // Acta Mathematica Scientia. 2009. V. 29. №6. P. 1629-1646.
12. Hirt C. W., Nichols B. D. Volume Of Fluid (VOF) method for the dynamics of free boundaries // J. Comp. Phys. 1981. V. 39. №1. P. 201-225.
13. Osher S., Sethian J. Front propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations // J. Comp. Phys. 1988. V. 78. №2. P. 12-49.
14. Yabe T., Xiao F., Utsumi T. The constrained interpolation profile method for multiphase analysis // J. Comput. Phys. 2001. V. 169. №2. P. 556-593.
15. Karni S. Multicomponent flow calculation by a consistent primitive algorithm // J. Comput. Phys. 1994. V. 112, №1. P. 31-43.
16. Abgrall R., Karni S. Computations of compressible multifluids // J. Comp. Phys. 2001. V. 169. №2. P. 594-623.
17. Yabe T., Wang P. Y. Unified numerical procedure for compressible and incompressible fluid // J. Phys. Soc. Japan. 1991. V. 60. №7. P. 2105-2108.
18. Fedkiw R., Liu X.-D., Osher S. A general technique for eliminating spurious oscillations in conservative schemes for multiphase and multispecies Euler equations // Int. J. Nonlinear Sci. Numer. Sim. 2002. V. 3. №2. P. 99-106.
19. Harlow F. H., Welch J. E. Numerical calculation of time-dependent viscous incompressible flow with free surface // Phys. Fluids. 1965. V. 8. P. 2182-2189.
20. Amsden A. A., Harlow F. H. A simplified MAC technique for incompressible fluid flow calculations // J. Comput. Phys. 1970. V. 6. №2. P. 322-325.
21. Patankar S. V., Spalding D. B. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows // J. Heat Mass Transfer. 1972. V. 15. P. 17871806.
22. Harlow F. H., Amsden A. A. Numerical simulation of almost incompressible flow // J. Comput. Phys. 1968. V. 3. №1. P. 80-93.
23. Issa R. I. Solution of the implicitly discretised fluid flow equations by operator-splitting // J. Comput. Phys. 1986. V. 62. №1. P. 40-65.
24. Ogata Y., Yabe T. Shock capturing with improved numerical viscosity in primitive Euler representation // Comput. Phys. Commun. 1999. V. 119. №2-3. P. 179-193.
25. Doihara R., Takahashi K. Numerical calculation of laser-produced bubble near a solid boundary until the second collapse // JSME. Int. J. Series B. 2001. V. 44. №2. P. 238-246.
26. Kawasaki K. Numerical model of 2-D multiphase flow with solid-liquid-gas interaction // Int. J. Offshore Polar Eng. 2005. V. 15. №3. P. 198-203.
27. Peng G., Ishizuka M., Hayama S. An improved CIP-CUP method for submerged water jet flow simulation // JSME Int. J. Series B. 2001. V. 44. №4. P. 497-504.
28. Ida M. An improved numerical solver for the unified solution of compressible and incompressible fluids involving free surfaces. II. Multi-Time-Step integration and applications // Comput. Phys. Commun. 2003. V. 150. №3. P. 300-322.
29. Nomura S., Nishida K. Numerical simulation of a single bubble rising in an ultrasonic standing wave field // Jpn. J. Appl. Phys. 2003. V. 42. P. 2975-2980.
30. Tong M., Browne D. J. Modelling compressible gas flow near the nozzle of a gas atomiser using a new unified model // Computers and Fluids. 2009. V. 38. №6. P. 1183-1190.
31. Tanaka N., Maseguchi R., Ogawara T. Improvement of conservative property and interface resolution of mesh-based two-phase flow simulation algorithms for splashing fluid behavior // Nuc. Eng. Design. 2010. V. 240. №12. P. 3984-3990.
32. Lee S.-J., Cho B.-G., Lee I. Two-dimensional unsteady aerodynamics analysis based on a multiphase perspective // Computers and Fluids. 2012. V. 53. P. 105-116.
33. Xiao F. A class of a single-cell high order semi-Lagrangian advection schemes // Month. Weath. Rev. 2000. V. 128. №4. P. 1165-1176.
34. Yabe T., Xiao F. Description of complex and sharp interface with fixed grids in incompressible and compressible fluid // Comp. Math. Applic. 1995. V. 29. № 1. P. 15-25.
35. Аганин А. А., Ильгамов М. А., Малахов В. Г., Халитова Т. Ф., Хисматуллина Н.А. Ударное воздействие кавитацион-ного пузырька на упругое тело // Ученые записки Казанского университета. Серия физико-математические науки. 2011. Т. 153. Кн. 1. C. 131-146.
36. Аганин А. А., Ильгамов М. А., Халитова Т. Ф. Ударное воздействие струи на жесткую стенку // Актуальные проблемы механики сплошной среды. К 20-летию ИММ КазНЦ РАН. T. 1. Казань: Фолиант, 2011. C.134-146.
37. Чижов А. В., Шмидт А. А. Высокоскоростной удар капли о преграду // ЖТФ. 2000. T. 70. Вып. 12. С. 18-27.
Поступила в редакцию 30.05.2013 г.
