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

Вычислительные технологии

Том 7, № 6, 2002

ОБ ОДНОМ СПОСОБЕ РАСЧЕТА БЕЗУДАРНОГО СИЛЬНОГО СЖАТИЯ ДВУМЕРНЫХ ГАЗОВЫХ

СЛОЕВ

С. П. Блутин, А. В. Рощупкин Уральский государственный университет путей сообщения,

Екатеринбург, Россия e-mail: SBautin@math.usart.ru

Numerical algorithm for two-dimensional shockless power compression is presented.

Results of some numerical experiments are described.

Математическое описание процесса безударного изэнтропического сжатия идеального газа до любого наперед заданного значения плотности, в том числе до бесконечной плотности (подробную библиографию см. в [1]) представляет интерес в связи с проблемой лазерного термоядерного синтеза. Особый интерес представляет сжатие неодномерных слоев газа, так как такое сжатие может оказаться энергетически более выгодным, чем одномерное как, например, (см. работы [2, 3]) при безударном сильном сжатии специальных призматических объемов. При согласовании угла призмы и показателя политропы 7 могут быть получены большие степени кумуляции газа, чем при сжатии одномерных объемов. К сожалению, такая конфигурация сжатия неустойчивая, так как малые изменения угла призмы от специального значения могут привести к возникновению ударных волн.

1. Постановка задачи

Рассматривается система уравнений газовой динамики в случае двумерных нестационарных изэнтропических течений политропного газа

(7 - 1)

аг + V ■ Уа + 1 [ ' а divV = 0,

Vt + гс^ х V + 1У (V ■ V) + Т-^ГгаУа = 0,

2 (7 - 1)

где а = р(т-1)/2, р — плотность, 7 — константа в уравнении состояния политропного газа р = р7/7, 7 > 1, р — давление; V = (м1,м2) — вектор скорости газа. Пусть в дальнейшем "О = (а,и1,и2) — вектор неизвестных функций.

Требуется найти решения системы (1), реализующие в момент времени £ = сжатие некоторого слоя газа до наперед заданной плотности р = р* (ж) в окрестности заданной кривой С*.

© С. П. Баутин, А. В. Рощупкин, 2002.

В системе (1) стандартным образом введены безразмерные переменные — за масштабные значения плотности, скорости и расстояния взяты соответственно плотность исходного однородного покоящегося газа, скорость звука c в этом газе и наименьший радиус кривизны заданной линии C*.

Для решения задачи о безударном сильном сжатии двумерных газовых слоев вместо декартовых переменных ж^ж2 вводится специальным образом выбранная ортогональная криволинейная система координат £,n. Если C* — кривая, в окрестности которой реализуется сильное сжатие, а Uo — априори заданное фоновое течение, например, однородный покой с р = р° = 1, которое сжимается до заданной плотности р*, то за координатные линии £ = const берутся прямые с направляющими векторами

V*(£) = V°(£) +

2

(Y -1)

a°n*(£), U0 = Uo|c _

где n*(£) — единичный вектор нормали к C*. Кривые, ортогональные прямым £ = const, есть координатные линии n = const.

Связь между декартовыми координатами x1,x2 и криволинейными ортогональными координатами £, n дается формулами

xi = Vi(£) - n , X2 = Ы£) + n , A(£) = vVi(£)]2 + (£)]2,

A(£)

X = ^(£), х2 = ^2(£) — параметрическое уравнение кривой п = 0. Якобиан такого преобразования имеет вид

J = -

1-

n

r(£ )J

A(£), r(£)

ii(£ V2(£) - 12(£М(£) A3(£)

-i

где г(£) — радиус кривизны кривой С* в точке £. Система (1) во введенных переменных имеет вид

Y - 1 Y - 1

at + + Ai(n,£)va + a К + Ai(n,£)vi] = A2(n, £)a,U

22 ut + uun + Ai(n'£)vuf + aan = -A2 ("'£)v'

2

Vt + uvn + Ai(n,£)vv? + Ai(n,£)-- aa? = A2(n,£)vu,

Y - 1

(2)

где

Ai (n,£)

r(£)

A(£)[r(£) - n]:

A 2 (n,£)

r(£) - n'

и, V — проекции вектора скорости на оси п и £ соответственно.

Согласно [1] решение задачи о безударном сильном сжатии получается при состыковке течений, являющихся решениями двух вспомогательных характеристических задач Коши (ХЗК). Первое течение — решение системы (2) с данными Коши на характеристике С( фонового течения и0

±

U(t, n, £) lc± = U°°(t, £) = U°(t, n, £) U

(3)

1

с одним дополнительным условием

да(£,п,С)

1

о, (4)

г=и

которое, во-первых, обеспечивает единственность решения первой ХЗК, а во-вторых, это условие придает решению первой ХЗК особенность, аналогичную особенности в центрированной волне Римана. Задача (2)-(4) в пространстве специальных переменных имеет единственное решение [1]. Второе течение является решением другой ХЗК с данными на характеристике С± первого найденного течения

И(£,п,С)1 с± = И^С) (5)

и дополнительным условием

а(£,п,С)|4=4, = а*(п,С) а1(£,С)|4=4, = а*(п,С)|п=п»(?) ' (6)

где С) — значение решения первой ХЗК на звуковой характеристике С±; п = п*(С) — кривая С* в координатах п,С В условии (6) а*(п,С) определяется значением р*(п,С) — наперед заданной плотности, до которой требуется сжать газ фонового течения. Задача (2), (5), (6) также имеет единственное аналитическое решение [4].

В случае аналитичности входных данных решение обеих ХЗК можно искать в виде специальных степенных рядов [1, 4]. В частности, для первой ХЗК при £ = £* решение дается формулами [1]

2 1 и(а, С) = и0 + ^ == ^а ± а0] ,

7 1 у 1 + [А1 (п* (С), С )п* ]2 ( С) о ± 2 А1(п*(С),С)п*

Па,С)= v ±-- . а,

7 1 у 1 + [А1 (п * (С), С )п* ]2 П1(а,С) = и Т\] 1 + [А1(п*(С), С)п *]2а - Л(п*(С),С)п*V,

где верхний знак берется при сжатии слоя газа изнутри, а нижний — при сжатии снаружи. Если найдены решения обеих ХЗК (2)-(4) и (2), (5), (6), то решение задачи Коши (ЗК)

Пр (£)|4=4о = П2 (С0) , ^СО^ = С0

дает в пространстве переменных £,п,С в параметрическом виде п = ПР(£),С = СР(£) траекторию частицы газа, которая в момент £ = £0 находилась в точке (п = п0(С0),С = С0), где п = п0(С) — некоторая кривая. Перебирая значения С0 из некоторого интервала [Е^, е2], мы получим в пространстве переменных (£, п,С) поверхность, которую можно считать поверхностью сжимающего непроницаемого поршня.

Области определения решений задач (2)-(4) и (2), (5), (6) представлены на рис. 1. Поверхность, ограниченная точками с номерами 9, 2, 3, 12, есть звуковая характеристика С± фонового течения; 5, 2, 3, 8 — звуковая характеристика С± решения первой ХЗК; 2, 3 —

кривая, на которой в момент £ = £* происходит сильное сжатие газа фонового течения, т. е. кривая С*. В момент £ = £* эта кривая есть область определения двумерного обобщения центрированной волны Римана; 1, 2, 3, 4 — область, в которой в момент £ = £* находится сжатый газ с заданным распределением плотности р*; 6, 1, 4, 7 — поверхность О1 — траектория движения сжимающего поршня по области течения, задаваемого решением второй ХЗК; 10, 6, 7, 11 — поверхность 02 — траектория движения сжимающего поршня по области течения, задаваемого решением первой ХЗК; 10, 1, 4, 11 — поверхность движения сжимающего поршня 01 + 02; область между поверхностями 9, 2, 3, 12 и 5, 2, 3, 8 — область определения обобщения центрированной волны; область между поверхностями 5, 2, 3, 8 и 1, 2, 3, 4 — область определения волны сжатия, непрерывно примыкающей к центрированной волне и задаваемой решением второй ХЗК.

Рис. 1. Области определения искомых течений газа.

Доказанные в [1, 4] теоремы утверждают, что существует ненулевая масса газа, которую можно безударно сжать до любой плотности. Однако эти теоремы не позволяют определить максимально возможную массу безударно сжимаемых двумерных слоев. Для определения этой массы можно применять численные методы решения двух поставленных ХЗК.

2. Численный метод

Будем считать, что фоновое течение U°(t, п,£) однородно и соответствует состоянию покоя. Тогда V°(£) = 0, а кривые C* и п = 0 совпадают.

Вместо u и а по аналогии с одномерным случаем [5, 6] введем новые неизвестные функции R и L:

2 2 R = u(t, п, £) +-тa(t, п, £), L = u(t, п, £)---a(t, п, £),

Y — 1 y — 1

а функцию v оставим без изменений. Рассмотрим кривые С+, С- и С°, лежащие в плоскости £ = const и удовлетворяющие уравнениям

-+ : d^ = Y + 1R _ Y — 3 L с- . dl = Y + 1 L _ Y — 3 R ' dt 4 4 , ' dt 4 4 ,

~о , ^ = Л + Ь ' ^ 2 '

Если бы искомое течение было одномерным (ди/д£ = 0), то введенные линии С+, С-, С° являлись бы соответствующими характеристиками системы (2) [5, 7].

Для функций Л, Ь и V система (2) после ряда эквивалентных преобразований примет

вид

Г dR

т

dL

и

dv dt

AivR? -

C+

Y-1

4

(R - L)vc - АЛ v2 - I—1 [R2 - L2]

AivL? (R - L)v? - AW v2 + Y 1 r°2 r2

C-

Y — 1

4

= -A1vvj +

Y - 1

8

Ai (R - L) (R? - ) + A2V

R2 - L2

R + L

(8)

с0

обозначает производную вдоль кривой C.

C

d

ГДе dt _

Для определенности реализацию численного метода будем рассматривать для случая сжатия газа изнутри и с конечным постоянным распределением плотности в момент сильного сжатия: а* = a(t *,£,n) = const.

Разобьем отрезок [51; 52] точками £k = 51 + кД£, Д£ = (S2 - 51)/N, и решение будем строить в каждой плоскости £ = £k, k = 0,1,..., N. Дополнительно предположим симметрию искомого течения при £ = £0, т. е.

д U

5£

= 0.

i=io

Для этого достаточно предположить симметрию линии С* относительно прямой £ = £0 и задать условие

да*

д£

0.

i=io

Рассмотрим построение разностной сетки и вычисление газодинамических параметров в ее узлах Рг для плоскости £ = £к, к = 0,1,N. Разобьем отрезок [1, а * ] на п частей

точками а0, г

(P0+1), г =

обозначения

0,

.П а°

а0, аП = а * так, чтобы угол между прямыми C (Pi) и

0,..., n - 1 был равен a/n, где а — угол между C (P°), и С (РП) и введем

с- (pfc): v = Gkt + и - Gktks ], G(t,n,£) = ^L(t,n,£) - ^R(t,n,£), Gk = G(Psk),

4

4

l = "г

Ln — "n -

Y - 1

а0г,

Pi = (M,£k) = Po = (t= 0,£k)

а0г - 1

R0 = u0 +

Y - 1

а0, u0

-2-

Y - 1

, v0 = 0, г = 0,..., n.

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

£k П C- = C7- (Po0) : V = G0t + [n0 - G0t0] .

8

2

2

Точки слоя, т.е. узлы сетки, будем обозначать Р3г = (£•, п3,, где г — номер слоя, г = 0,1,..., ] — номер точки на слое, ] = 0,1,...

Для построения новой расчетной точки Р. необходимо знать две точки Р3г—1 — предыдущую рассчитанную точку со слоя I и точку Р,-^ — точку, полученную при пересечении прямой С0 (Р3г—1), выходящей из точки Р._1, и I — 1-го слоя

С0 (р?): п = ф?£ + [п? — Ф?£к], ф(£,е,п) = Р(£,е,п) +ь(£,е,п), Ф? = Ф(р*).

Если Р]_1 = Р0, то р,_1/2 — первая точка слоя, выходящего из Р0-1, I > 0.

Получены следующие формулы для расчета узлов и параметров газа в них: координаты новой расчетной точки Р- находятся как координаты точки пересечения прямой С_ (Р3г—1),

выходящей из точки Р]_1, и прямой С+ ^Р,!_1/2^, выходящей из точки Ргг—1/2:

(Р?) : п = Р?£ + [п? — Р?£?] ,

р(£, п, е) = ^д(£, п, е) —^ь(£, п, е), р? = р (р?).

Тогда в первом приближении

£; = /-—1/2 о3_1 ; = ¿3_1/-—1/2 Рг_1/2о3_1 ¿3_1 Р-—1/2 ¿3_1 Р-—1/2

/г—1/2 = ''г_1/2 Рг_ 1/21/2, У?_1 = ¿3_1£3_1-

Значения функций Р и Ь в найденной точке Р. в первом приближении считаются по формулам

р3=рг—1/2+(£.—£-—!/2)я1 (р;-^) , 4=4_1+(£.—£. _ 1)^2 (р?_О ,

где Н1(Р), Н2(Р) — значения правых частей первого и второго уравнений системы (8) соответственно в точке Р. Через точку Р. проводим прямую С° (р) до пересечения со слоем I — 1, и пусть Рг;-1/2 — точка пересечения. Значения функций Р, Ь и V в точке Рг;—1 /2 определяются интерполяцией по точкам отрезка, содержащего нашу точку. Тогда значение функции V в точке Р. в первом приближении определяются следующими соотношениями:

V = гй/2 + (£. — 4-1 /2)Нз(рГ—1/2) ,

где Н1(Р) — значение правой части третьего уравнения системы в точке Р. Пересчет точки производится по формулам

71 _ 1 — я, 7 _ 1 — _ 1 71

£, = / г _ 1/2 73 _ 1 ; = ¿3 _ 1 /г _ 1/2 Рг _ 1/273 _ 1 3 = — Р _ 1 , п = ¿г — Р _ 1 ,

Рг_ 1/2 ¿3—1 Рг

т3 _1 ^ г_ 1/2 ^3 _1 ^ г_ 1/2

7; _ 1 = пг _ 1 — 7; _ 1 £г _ 1 7г = пг — ¿7; £г

/г_ 1/2 = пг _ 1/2 Рг _ 1/2£г _ 1/2, 73_ 1 = п3_ 1 ¿3_1£3_1,

р _ 1 = Р_1/2 + Р3 = ¿3 _ 1 + ¿3

г—1/2 2 ' 3 _1 2

Для пересчета значений функций Л, Ь и V нужно знать значения их производных функций по координате £ в точке Р?. Для этого на плоскости £ = 1, к > 0, нужно найти точку с координатами (Р?.¿, Р?.п,£к-1), где через Р.£, Р.п обозначены координаты £ и п точки Р. Значения искомых функций в этой точке нужно интерполировать по вершинам треугольника, содержащего данную точку. В предположении, что при достаточно малых Д£ = £к — £к-1 координаты узлов сетки на плоскостях £к и £к-1 отличаются не сильно, этот треугольник надо искать в окрестности узла плоскости £ = £к-1 со слоя I и номером ]. Тогда

^ R(P■n,£k-Q - Rj _ г ^ L(P ^Pj ■n,£k-Q - Lj

j) - Д£ , ?( j) - Д£ ;

v(jtj,£k-1) - vj v? (Pj) - ——-—2

и значения функций R и L пересчитываются по формулам

R = R-1 + (i - ¿1-1 ) u г-1/2 Rj = Ri-1 /2 + (tj ¿г—1/2)

H1(P_-i/2) + H1(Pj)

Lj = Lj-1 + (tj - tj-1)"

Н2(Р?_1) + Я2(Р]) 2 '

Пересчет точки Р,-^ производится следующим образом: находятся координаты точки пересечения прямой

ф'-1 + ф?

/о ^

2

,1-1

V = Ф1 + (nj - фjj , ф

r-1/2 ~ ^ j j > 2

выходящей из точки Р? и слоя I — 1. Пересчет значения функции V ведется по формуле

,.,' = +(г- ) Яз(Р-'-1/') + )

^ = ^-1/2 + V-1/2) 2 '

Я о'-1 _ о'-1

^-1 ,

Далее полагаем Pj-1 = Pj, Pi-/2 = Pr-1/2 и продолжаем расчет до тех пор, пока

возможно нахождение пересечения прямой С° (Pj) с предыдущим слоем l - 1.

В ходе расчета нужно следить за размерами получающегося треугольника с вершинами в точках Pj-1/2, Pj, P,-1/2. Если его стороны превосходят заранее определенную величину, то для сохранения точности расчетов вводится новый слой с начальной точкой, лежащей посередине между точкой Pj- Ц/2 и точкой Pj.

После построения разностной сетки в плоскости £ = £k определяется траектория точки, лежащей на поршне. Для этого надо рассмотреть разностное приближение задачи Коши (7). Начальное положение этой точки можно выбирать двумя способами: если мы хотим, чтобы в момент времени t = t* сжатый газ находился между линиями п = П2(£) : |П2(£) - П*(£)1 = const и п = П*(£), то начальное положение поршня надо задать в точке (¿*, П2(£k), £k), k = 0,..., N. При этом траектория поршня строится в сторону убывания времени и может получиться, что для различных точек на поршне построенная траектория придет на звуковую характеристику C- фонового течения в разные моменты времени.

Второй способ: начальное положение можно задать на звуковой характеристике C-фонового течения (t0,C*(t0,£k),£k), k = 0, ...,N. При этом траектория поршня строится в

2

сторону возрастания времени и может получиться, что |п2(£) _ п*(£)| = const — ширина итогового сжатого слоя не будет одинакова при разных значениях £.

После построения поршня производится вычисление массы несжатого газа m° и массы сжатого газа m*. При сравнении этих двух величин определяется точность расчета.

3. Результаты численных экспериментов

Для проверки правильности работы алгоритма проведен ряд тестовых расчетов, результаты которых сравнивались с одномерными расчетами по сжатию цилиндрически симметричных слоев газа [6]. Для этого, как и в одномерном случае, в качестве кривой C* брался круг единичного радиуса, а все производные по координате £ полагались равными нулю. Оба алгоритма на входе имели одинаковые параметры, в том числе а° = 1, а* = 1.5, At = 0.01. Результаты расчетов приведены в табл 1.

Здесь первая цифра в номере варианта указывает на расчеты по одномерной или двумерной программам, а вторая цифра соответствует следующим значениям: 1 — y = 1.4, n = 140; 2 — y = 1.6, n = 80; 3 — y = 1.7, n = 80. Из таблицы видно, что при малых а* результаты расчетов по одномерной и двумерной программам близки.

Также были проведены сравнения с квазиодномерными расчетами: в качестве кривой C* был взят эллипс с большой полуосью, равной 2, и малой полуосью, равной 1, а производные по координате £ на каждой плоскости £ = const полагались равными нулю. Один из результатов приведен в табл. 2.

Здесь вариант 4 — квазиодномерный случай, вариант 5 — расчеты в двумерном случае. В обоих вариантах а° = 1, а* = 2, y = 1.7, N = 30, n = 150, At = 0.01. В варианте 5 — а° = 1, а* = 2, y = 1.7, что соответствует р° = 1, р* = 6.89. Заметим, что бесконечно сильная ударная волна для данного показателя y дает скачок плотности от р° = 1 до р* = 3.87 [7]. Газодинамические параметры между поршнем и фоновым течением монотонны как по пространственной, так и по временной координатам. Максимальные значения р и |V| достигаются на поршне. В варианте 5 максимальное значение |V| = |(u, v)| = \J(_ 2.92719)2 + (0.00097)2 = 2.92719, максимальное время движения поршня 0.1.

Таблица 1 Результаты сравнения с одномерными расчетами

№варианта m° m* Am (Am/m°)100

1.1 1.7731148937 1.8015651198 0.0284502261 1.6045337051

2.1 1.7731148937 1.7955082644 0.0223933707 1.2629396312

1.2 1.9726060272 1.9983039305 0.0256979033 1.3027387601

2.2 1.9726060272 1.9818015414 0.0091955142 0.4661607069

1.3 2.1208892004 2.1417917549 0.0209025545 0.9855561760

2.3 2.1208892004 2.1275818075 0.0066926071 0.3155566589

Т а б л и ц а 2 Результаты сравнения с квазиодномерными расчетами

№варианта m° m* Am (Am/m°)100

4 0.9398832612 0.9811644932 0.041281232 4.3921658676

5 0.9398832612 0.9811657192 0.041282458 4.3922963444

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

Расчетами установлено, что в тех плоскостях £ = £к, где кривизна линии С* на каких-либо промежутках изменяется не резко, расчеты двумерного течения достаточно хорошо совпадают с квазиодномерными расчетами. Данный факт, скорее всего, является следствием соответствующего выбора криволинейной системы координат, связанной с линией С*. Основная погрешность вычисления в обоих случаях была сосредоточена в плоскостях, где изменение радиуса кривизны происходит быстро.

Четверть поверхности движущегося поршня для варианта 5 графически представлена на рис. 2.

Расчеты течений до больших степеней сжатия требуют существенного уменьшения размеров расчетных ячеек. Например, при расчете уже одномерных течений в работе [6] число характеристик, выходящих из особой точки, порядка 100 000. Уменьшение размеров расчетных ячеек приводит к существенному увеличению времени счета.

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

Список литературы

[1] БАУТИН С. П. Математическая теория безударного сильного сжатия идеального газа. Новосибирск: Наука, 1997.

[2] Сидоров А. Ф. Некоторые оценки степени кумуляции энергии при плоском и пространственном сжатии газа // Докл. АН СССР. 1991. Т. 318, №3. С. 548-552.

[3] Сидоров А. Ф. Оценки предельных степеней кумуляции энергии при безударном сжатии газа // Докл. АН СССР. 1993. Т. 329, №4. С. 444-448.

[4] БАУТИН П. С. Двумерные нестационарные задачи о получении наперед заданных распределений параметров течения идеального газа. Екатеринбург: УрГУПС, 2000. Деп. в ВИНИТИ от 10.04.2000 за №968-В00.

[5] баутин с. п., Николаев ю. в. Об одном методе расчета безударного сильного сжатия одномерных слоев газа // Вычисл. технологии. 2000. Т. 5, №4. С. 3-12.

[6] николаев ю. в. О численном решении задачи безударного сильного сжатия одномерных слоев газа // Вычисл. технологии. 2001. Т. 6, №2.

[7] рождественский б. л., Яненко н. н. Системы квазилинейных уравнений и их приложения к газовой динамики. М.: Наука, 1978.

Поступила в редакцию 11 сентября 2001 г., в переработанном виде — 9 июля 2002 г.