CC BY
103
НЕЯВНАЯ СХЕМА С.К. ГОДУНОВА ПОВЫШЕННОЙ ТОЧНОСТИ ДЛЯ РАСЧЕТА СТАЦИОНАРНЫХ СВЕРХЗВУКОВЫХ ТЕЧЕНИЙ
М. Я. Иванов, Р. 3. Нигматуллин, Л. В. Терентьева
В последние годы предложено несколько разностных схем для интегрирования уравнений Эйлера, позволяющих рассчитывать течения со скачками уплотнения методом сквозного счета. Библиография по этому вопросу имеется, например, в статье /2/. Особый интерес представляют те из схем, которые обеспечивают описание скачков уплотнения и других разрывов газодинамических параметров при отсутствии паразитических колебаний решения. Этот интерес обусловлен тем фактом, что в условиях наблюдаемого в последнее время быстрого прогресса вычислительной техники и методов автоматического построения сеток (в том числе и адаптивных) недостаток надежных и точных разностных схем является главным препятствием дальнейшему развитию вычислительной аэрогидродинамики.
В областях со сверхзвуковой составляющей скорости вдоль направления, выбранного в качестве маршевого, использование "бесконечно большого" шага по времени позволяет в некоторых случаях свести интегрирование уравнений Эйлера к решению системы уравнений, описывающей стационарные движения /3/. В настоящей работе на основе неявного варианта схемы Годунова повышенной точности, разработанного в /2/, построена неявная схема для расчета стационарных сверхзвуковых течений маршевым методом. С целью повышения точности аппроксимации производных по поперечной координате используются кусочнопараболические распределения параметров по ячейкам разностной сетки, удовлетворяющие условиям монотонности. Расчет потоков на границах ячеек проводится с помощью процедуры распада произвольного разрыва, учитывающей локальную струюуру потока. Обращение неявного оператора проводится также с учетом местной картины течения посредством анализа собственных значений матриц Якоби.
В качестве примеров приведены расчеты течения невязкого идеального газа в плоском канале и истечения осесимметричной недорасширснной струи в затопленное пространство. Проведено сопоставление полученных результатов с опубликованными расчетными данными.
1. Разностная схема
Система уравнений Эйлера, описывающая стационарные плоские (V =0) и осесимметричные ( у=1) течения невязкого идеального газа, может быть записана в консервативной форме в декартовой (или цилиндрической) системе координат (х,у) следующим образом:
ри " V ’ 0 ‘
где Р = р + ри2 , 0 = риу , f= 0
риу р + ру ур
_(е + р)и_ (е + р)у 0
Здесь х, у - декартовые (или цилиндрические) координаты, р - давление, р - плотность, и, V - компоненты скорости вдоль направлений х, у соответственно, е - полная энергия единицы объёма газа, связанная с газодинамическими переменными уравнением состояния:
Все переменные представлены в безразмерном виде. Скорость отнесена к критической скорости Укр , плотность - к критической плотности ркр , давление - к
удвоенному критическому скоростному напору - ркрУ2р, координаты - к
линейному размеру задачи Ь.
_ » дБ _ 90 т.
Определим матрицы Якоби А =—, В = —, где II - вектор
эи
независимых переменных и = (р, и, V, р) .
Для "х-сверхзвуковых" течений система (1.1) является гиперболической, матрица А - невырожденная, и система (1.1) может быть записана в виде
а^+А-'в^А-’г.
Эх Эу
Матрица А-1В имеет четыре собственных значения >М 2 = 1&0. 'Ч = 1ц(0 + ц), Я,4 = 1д(в - ц), где 1§0 = у/и, с1£ц = л/м2 -1, М - число Маха, и может быть представлена в виде А~!В = БАЗ-1, где Л = сНа§(Х.)Д2''Ч'^4) ■
диагональная матрица, Б - матрица преобразования, приводящего А-1В к диагональному виду.
При построении неявной процедуры интегрирования по маршевой координате х будем исходить из следующего соотношения:
в . 1
5Рп-! +•
. ар"
-1г--------+ -
г , эр —1г
п-1
1 -ь в + г дк 1 + в + г
Эх
+о|(1 +• 2я+Зг - 2q)h2 +(1-2г-Зч)Ь3 +Ь4|.
(1.2)
где ЭРП =Рп+1-Рп .
Верхний индекс отвечает своему значению координаты л\ х=п11, где Ь -шаг интегрирования по х. Формула (1.2) задает двухпарамстрическое семейство схем второго порядка по .V (при выполнении равенства 1 + 2з + Зг-^=0) и однопараметрическое семейство схем третьего порядка при выполнении равенств
Ч = (1-2г)/3, 8 = -(1+13г)/б. (1.3)
Равенство (1.2) дает следующую неявную схему для интегрирования системы (1.1). Для простоты рассмотрим плоский случай:
1 + в + г ду 1 + б + г 1+Б + Г ду И-Б + Г ду
Для оценки величины-------- воспользуемся тем, что 50 п =ВА 5РП, а
ду
матрица ВА-1 может быть представлена в виде ВА-1 = а(а-1в)а-1 = А8А8-1А-1 = ВАБ-1, где 55 = АЭ .
Выполним расщепление матрицы ВА-1 на две составляющие с учетом локальной характеристической картины течения. Представим ВА-1 в виде разности ВА-1 = (}+ - СГ Здесь
А± = (|А| + л)/2,
|л| = Шая(|^|,|?41ЫМ)-
Обе введенные матрицы С>+ и 0" имеют неотрицательные собственные значения. Неявная консервативная разностная схема имеет вид
----ВА-11бРп =ЯП,
1 + Б + Г ду _|
где И11 - известный на п шаге по х вектор правых частей
Кп=__1_5РпЧ---------------------------
1 + в + Г 1 + Б + Г ду 1 + Б + г ду
Переход от дифференциального оператора в левой части осуществляется с помощью односторонних разностей, записываемых с учетом знаков собственных значений. Имеем
1+—3—ь(уд+-до~)
1+8+Г V ' /
8РП=Я\ (1.4)
где разностные операторы определены формулами УГj = ^
Правая часть (1.4) рассчитывается по параметрам с предыдущих слоев по координате х. Для оценки производных по у использовалась процедура распада разрыва С.К. Годунова [1]. Для повышения порядка аппроксимации задавалось кусочно-параболическое распределение параметров в ячейках, удовлетворяющее условиям монотонности. Подробное изложение метода дано в [4] на примере скалярного одномерного нестационарного уравнения. В случае стационарной системы уравнений параболическое распределение по ячейкам в направлении у строится для "характеристических” переменных определяемых посредством соотношения 8\У = 8-18и.
Отметим, что при выполнении равенств (1.3) схема (1.4) имеет третий порядок лишь в том случае, когда матрицы С>+ и СГ вычисляются на слое п +1/2
(если их вычислять на п-том слое, порядок схемы понижается до второго). Для того чтобы устранить этот недостаток, можно, например, на каждом шаге интегрирования по х ввести дополнительные итерации по схеме
1+—3—ь(уд+-дд_)15Р^ = кН (1.5)
. _ 1 + Б + Г V
§р(р) _ р(р+1) _р(р); р(°)_рп^ рп+1_р(р+1)
я(р)=Ё:п+в(р), ёп=к"+—а_
1 + Б + Г ду
0(р) = -5р(рЧ)-----!—ь в^ = 0.
. 1 + Б + г ду
Если р=1, то схема (1.5) совпадает со схемой (1.4).
Для решения системы уравнений (1.4) необходимо обращение блочных трехдиагональных матриц, что реализовано методом матричной прогонки.
При решении многих прикладных задач удобно вместо декартовых координат X, у использовать криволинейные координаты. Перепишем систем)' (1.1) в криволинейных координатах £ = ^(х,у), г| = Л^у) с сохранением дивергентной формы уравнений:
* (1.6)
где Р = 1_1(р§х + С^у)у \ 6 = + СЛу)у \ ! = Г*Г.
Здесь I - якобиан преобразования координат 1 = £хг|ч. - ^уг|у.
Система (1.6) может быть записана в виде
8е, дц \ * I
* др л ай
Матрицы Якоби А = —, В = — связаны с матрицами А и В
аи эи
соотношениями А = 1(А£Х + В£у), В = 1(Аг)х + Вг|у), а собственные числа
Л | А
матрицы А В определяются следующим образом
^ 1 = 1’2’3>4- 0^
Неявная консервативная разностная схема для интегрирования системы (1.6) строится аналогично изложенному выше. Разностные уравнения для ]-й расчетной ячейки имеют вид:
1+_Л—ь(у<3+-ЛСг) =£/"), (1.8)
. 1 + 8 + Г '
<3* = ЗА* Г1, А* = (|Л| ± А)!%
В настоящей работе для всех вариантов построенной неявной схемы использовался неявный способ задания граничных условий с использованием информации о локальном угле наклона характеристик. Все необходимые сведения о постановке граничных условий в неявном виде можно найти в 15/.
2. Результаты расчетов
Приведем некоторые численные результаты, демонстрирующие работоспособность различных вариантов разностной схемы (1.8). Первый результат относится к истечению недораспшренной осесимметричной струи в затопленное пространство при п=р0/ре=2 и М0=3, где ре - давление в затопленном
пространстве, а индекс нуль приписан параметрам на срезе сопла.
При численном решении этой задачи (как и следующей) удобно использовать следующее преобразование координат
1--7Т#? Ш)
у+(х)-у_(х)
где У-(х), у+(х) - уравнения, описывающие форму верхней и нижней границ
рассчитываемой области. В данном случае нижней границей является ось симметрии, а верхняя граница определяется в процессе расчета по заданному давлению в окружающем пространстве. Во входном сечении на срезе сопла задавался однородный поток, параллельный оси х. Число ячеек в поперечном направлении (NJ) равно 40. Параметры схемы: q = 0, s = 0, г = 0. Поскольку схема
является явной, то шаг по х ограничен условием устойчивости
1.
—тахЫ £1, i = 3,2,3,4, j = l,...,NJ.
Ду i,j 1 1
Расчеты проводились при числе Курангга, равном 0,2.
На рис. 1 приведены изобары поля течения. Здесь масштабы по осям х и у разные, а линии р = const даны через Др = 0,002. На рис. 3 нанесены кривые распределения давления в поперечных сечениях струи (у0 - отношение ординаты к ординате границы струи, координата х отнесена к радиусу сопла). Пунктиром нанесены распределения давления в тех же сечениях из работы [1], в которой расчет проводился явным методом Годунова первого порядка точности по х и по>\ Отметим, что схема с кусочно-параболическими распределениями параметров в ячейках дает более четкие профили скачков без осцилляций даже при первом порядке по маршевой координате.
Следующий результат относится к течению газа в плоском сужающемся канале. Угол наклона нижней границы равен 10° , после того как ширина канала составит 2/3 от ширины на входе, нижняя граница становится горизонтальной. На входе задавался равномерный поток, параллельный оси х, число Маха набегающего потока равно 3. Использовалось преобразование координат (2.1). Проведены расчеты с различными параметрами схемы. На рис. 2 приведено поле изобар течения, рассчитаное при q = 1/3, s = -1/6, г = 0 и числе Куранта, равном 1.0. Число дополнительных итераций на каждом шаге интегрирования было равно 3.
Линии р = const даны через Др = 0.015. На рис. 4 показано распределение давления вдоль нижней стенки канала. Пунктиром показаны результаты расчета по схеме первого порядка точности с q = 0, s = 0, г = 0 и числом Куранта 0.2. Сравнение результатов позволяет сделать вывод, что оба метода обеспечивают высокую точность расчета, однако использование неявных методов делает возможным применение больших шагов по х.
Рис.1
Рис.2
ооо гоо .чоо боо аоо у* і ооо Рис. 3
' 50
Рис. 4
СПИСОК ЛИТЕРАТУРЫ
1. Годунов С.К., Забродин А.В., Иванов М.Я. и др. Численное решение многомерных задач газовой динамики. М.: Наука, 1976.
2. Иванов М.Я, Нигматуллин Р.З. Неявная схема С.К. Годунова повышенной
точности для численного интегрирования уравнений Эйлера // ЖВМиВФ. 1987. Т. 27, № 11. С. 1725-1736. ^
3. Чакраварти С.Р., Жем К.-Й. Расчет трехмерных сверхзвуковых течений с дозвуковыми зонами на основе уравнений Эйлера // Аэрокосмическая техника. 1987. № 11. С.22-34.
4. Иванов М.Я., Крупа В.Г., Нигматуллин Р.З. Неявная схема С.К. Годунова повышенной точности для интегрирования уравнений Навье-Стокса // ЖВМиВФ. 1989. Т. 29, № 6. С.888-901.
5. Чакраварти С.Р. Уравнения Эйлера - неявные схемы и граничные условия // Аэрокосмическая техника. 1984. № 2. С.58-67.
ГИДРОДИНАМИЧЕСКИЕ ХАРАКТЕРИСТИКИ ВЫСОКОВОЛЬТНОГО
ЭЛЕКТРИЧЕСКОГО РАЗРЯДА В ДВУХФАЗНЫХ СРЕДАХ
В.Г.Ковалев
При обработке материалов электрогидравлическим способом в качестве рабочей среды, в которой развивается канал разряда и посредством которой передается возмущение, используется вода, которая, за исключением случаев специального эксперимента, содержит пузырьки газа. Содержащиеся в воде газовые включения, практически не изменяя плотности несущей среды, могут значительно изменить ее сжимаемость, в результате чего могут существенно измениться условия развития плазменного канала высоковольтного электрического разряда в воде и, как следствие, его гидродинамические характеристики.
Имеются немногочисленные работы, посвященные исследованию подводного взрыва ВВ в пузырьковых жидкостях, которые показывают, что присутствие в жидкости пузырьков газа оказывает заметное влияние на излучаемую при взрыве волну давления.
Однако взрыв ВВ - явление более простое, чем высоковольтный разряд в жидкости. Ударная волна полностью формируется в теле заряда ВВ независимо от характеристик окружающей среды. Напротив, при подводном электроразряде весь процесс выделения энергии существенно зависит от механических свойств окружающей среды, которыми определяется давление плазмы в канале разряда и тем самым - удельная электропроводность плазмы. Давлением плазмы, т,е. тем
