CC BY
38
Вычислительные технологии
Том 10, № 2, 2005
ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ АВТОКОЛЕБАТЕЛЬНЫХ ТЕЧЕНИЙ НЕЯВНОЙ СХЕМОЙ ЧЕТВЕРТОГО ПОРЯДКА*
В. И. Пинчуков
Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: [email protected]
Unsteady flows of an ideal gas with shocks and contact discontinuities are investigated. Two-dimensional Euler equations are solved by implicit fourth order Runge — Kutta scheme. The results of numerical calculations are presented.
Введение
Развитие численных методов и рост мощности ЭВМ позволяет использовать методы математического моделирования для изучения все более сложных течений в аэродинамике. Особый интерес вызывают автоколебательные течения ввиду их недостаточной изученности и опасности нестационарных режимов для аэродинамических конструкций. Автоколебательные течения изучались экспериментально и математически в [1-7]. По нашему мнению, расчет автоколебательных течений предъявляет особые требования к численным методам.
Чтобы точно отделять стационарные и нестационарные режимы, схемы для расчета таких течений должны сочетать противоречивые свойства, а именно, они должны быть малодиссипативными и не подавлять нестационарность малой интенсивности вблизи границы раздела режимов, тем не менее они должны обеспечивать надежную сходимость по времени в случае стационарных режимов сложной структуры. Кроме того, зачастую динамика автоколебательных течений определяется наличием в них отрывных зон. Учитывая, что размеры и структура этих зон весьма чувствительны к величине схемной вязкости, целесообразно использовать методы высоких порядков. Наконец, ввиду разрывного, как правило, характера автоколебательных течений также необходимо хорошее разрешение разрывов.
Всеми указанными свойствами обладают две неявные консервативные схемы Рунге — Кутты четвертого порядка [8, 9] (одна со скалярными, другая с матричными стабилизирующими операторами). Пространственный шаблон этих схем (в том числе операторы искусственной диффузии) построен с использованием элементов адаптации, что позволило улучшить разрешение разрывов и сходимость по времени. В данной работе предпринята
""Исследования, выполненные В.И. Пинчуковым и описанные в этой статье, стали возможными отчасти благодаря поддержке (грант RM1-2324-NO-Ü2) Американского фонда гражданских исследований и развития для независимых государств, образованных на территории бывшего Советского Союза (CRDF).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2ÜÜ5.
попытка исследования автоколебательных течений на основе матричной схемы, оказавшейся, согласно результатам тестовых расчетов [9], более точной при числах Куранта порядка единицы и выше.
1. Основные уравнения и численный метод
Двумерные осесимметричные уравнения совершенного невязкого газа могут быть записаны в дивергентной форме двумя способами. Для каждой из них была построена неявная матричная схема Рунге — Кутты четвертого порядка. Более предпочтительной в расчетах оказалась следующая запись:
(уи )4 + (уР )х + (уС)у = Н, (1)
U = (р, e, pu, pv), F(U) = (pu, (e + p)u, pu2 + p, puv), G(U) = (pv, (e + P)v, puv, pv2 + P), H(U) = (0, 0, 0, P),
e = p[e + (u2 + v2)/2], e(p, p) = P/[(k — l)p],
где p — плотность; u, v — компоненты скорости вдоль осей x и y; P — давление; e — полная энергия; е — удельная внутренняя энергия. Показатель адиабаты во всех задачах принят равным к = 1.4. Поскольку используются криволинейные сетки, произведен переход к произвольным неортогональным координатам, в результате чего система принимает вид
(yIU)t + (yFyv — yGxv+ (yGx? — yFy?)n = H1, I = yvx? - y?xn.
Уравнения численно интегрируются в безразмерных переменных: p = p p^, P = P P^, u = u (Prx/prx)l/2, v = v (P'x/p^)1/2, x = x l, y = y l, где l — характерный размер в задачах.
Приводимая ниже схема отличается от описанной в [9] лишь наличием дополнительных членов, аппроксимирующих правую часть уравнения (1). Дискретизация по времени правой части производится так же, как в неявной схеме Рунге — Кутты четвертого порядка [10] для обыкновенного уравнения, абсолютная устойчивость которой доказана с помощью линейного анализа. Согласно [10] в стабилизирующий оператор добавляются члены, которые в результате процедуры факторизации приводят к появлению новых множителей — стабилизирующих операторов. Факторизация позволяет свести обращение пятиточечного оператора с матричными коэффициентами к совокупности скалярных прогонок и локальных матричных обращений. Следует отметить, что описываемые схемы реализованы для случая произвольных криволинейных неортогональных координат. Однако, чтобы не усложнять изложение, здесь рассматриваются уравнения (1) в декартовых координатах. Используются обозначения 5±u = ±[u(x1, ...,xi ± Axl,...,xm) — u(x1, ...,xm)], A±u = 8± u/Axl.
Итак, рассмотрим четырехшаговую схему
yik
' r DH frl\4
1 — 4yDU + AUJ
т DF т
y¡k + yjk 4 DU (Sx )-1Лх Sx — -(sx)-1V:Sx
1
х
• т DF т
vL + vL 4 ^ (s; )-1Лу s ; - ; )-1v;s ;
u(1) — Un
т/2
I Г)П TTU _ T ГПТТП ,
+ B:k — H:k — V u:k ;
(2)
7 v:k
т DH {т!\
1 — 4VDU + HiJ
• т DF т
vL + vL 4 DU (Sx )-1Лх sx — 2(sx)-1v;sx
х
х
vL + vL
т DF т
3 4m(sy)-1луs; — 2(s;)-1v;s;
TT(2) — Г/(1) TT(1) — Г/ (n)
Uik Uik + v Uifc Uifc +
т/2
т/2
+ (Bk + B(1) — Hk — H^ )/2 — Vn(Uik + U^ )/2,
(1)
(1)
(3)
7
vik
. т1
1 + n -v
TT(3) tTn
vL + ет4Д4х — 2v;) (vL + £т4я% — 2v;) ik — ik+
+B® — Hk
V nUk,
(4)
7
vik
, т1
1 + n
v
vLk + É^R* "V
т
vik + — 2 v;
US+1 — u
ik
+ (B:nk + 4B((2) + — Hk — 4H^ — HT )/6 — Vn U:
>(3)
(2)
(3)
nn 4k •
(5)
Здесь — матрицы, приводящие к диагональному виду якобианы DF/DU, DG/DU:
Sx(DF/DU)(Sx)-1 = — диагональная матрица; l — максимальное собственное
значение матрицы DH/DU (в данном случае имеется лишь одно собственное значение Л = v(k — 1)); Лх, Лу — операторы дифференцирования четвертого порядка:
лх = д2(1 — №/6), лу = Д0(1 — №/6);
В имеет смысл баланса потоков для элементарной разностной ячейки; Д4ж, Д4у — вспомогательные операторы (см. далее); V = МХ + Уу, МХ, Му — нелинейные адаптивные операторы искусственной диффузии, переключающиеся в зависимости от гладкости решения (диффузия пятого порядка на "гладких" решениях и первого возле ударных волн):
^ = А-[(У' е)<+> + е^+к ^,
4
1
4
1
4
1
= A-[(yjе)^ - 4(yjе(6;))i+1/2fc№№, e(6x) = ax72b, e(2x) = ax(1 - 7)166. (6)
Здесь ax = c + |u|, c — скорость звука; b = 1/32, V;.*, Vy* — "урезанные" пятиточечные диффузионные операторы, используемые в стабилизирующих операторах; 0 < Ytfc < 1 — индикаторы локальной гладкости, вычисляются по давлению на основе формул
Yifc = 1 - [1 - max(Yik, 0)]2, Y°fc = [min(agifc, 1) - c]/(1 - c), (7)
дгк = (г-г+ + ерк)/{шах[(г+)2, (г-)2] + 10-9}, г- = 6-Ргк, = 6+Ргк.
Константы а и с позволяют регулировать "чувствительность" величин ^¿к к значениям разностной функции Ргк в узлах сетки, параметр е предназначен для предотвращения понижения порядка диффузии на гладких экстремумах или участках почти постоянного решения. По результатам пробных расчетов были выбраны следующие значения эмпирических величин: с = 0.4, а =1.3, е = 10-3.
Адаптивная диффузия (6) переключается в зависимости от локальных индикаторов "гладкости", она увеличивает диссипацию возле ударных волн. Формулы (7) обусловливают для слабых разрывов квадратичную зависимость коэффициента вязкости низкого порядка от градиента давления. Это обеспечивает хорошее разрешение слабых разрывов.
Вспомогательные операторы Д4х, Я4у также строятся по адаптивным формулам:
П4х = Д-[Д+(аХк)4ук7гкА- - (ах+1/2к)4У+1/2к(1 - Т*+1/2к)4/Дх2]Д+.
В соответствии с потребным порядком аппроксимации схемы для расчета величины В, имеющей смысл баланса потоков для элементарной разностной ячейки, можно использовать формулу четвертого порядка В = Лх(у^)гк + Лу(уС)гк. Однако при расчете течений с сильными ударными волнами вместо этой формулы целесообразно применение формул с элементами адаптации. По результатам пробных расчетов наилучшими оказались формулы
Вп = Д-^+1Ак + А-ёПк+1/2, ^+1/2,к = [^Хк + (^)П+1,к]/2 - Тг+1/2,к(/¿+1,к + /¿,к)/12,
Дк = шШ(7г+1/2,к,7г-1/2,к^к, 7г+1/2,к = шт(7г+1,к, 7г,к),
которые эквивалентны исходной формуле четвертого порядка, если индикаторы гладкости 7 принимают единичные значения. Их применение позволяет избежать появления отрицательных значений давления и плотности и других аварийных ситуаций, а также улучшить разрешение разрывов.
Согласно линейному анализу [10], дискретизации по времени схемы (2)-(5) дифференциальных и источниковых членов уравнений (1) абсолютно устойчивы (для многомерного уравнения переноса в случае нефакторизованной схемы и обыкновенного уравнения первого порядка), если
С = (т250/243)3/384 = (2 ■ 250/243)3/384, Л =1/(3 ■ 163), п =1.4106 ■ 10-3
(т — число пространственных переменных). Анализ совместного использования дифференциальных и источниковых дискретизаций в одном уравнении не производился, однако проблем с устойчивостью для схемы (2)-(5) в расчетах не наблюдается, если выполнено эмпирическое ограничение К < 3 (К — число Куранта). Напомним, эта же схема для случая плоских течений работоспособна [9, 10] при ограничении К < 5.5, которое возникает, по-видимому, вследствие нелинейного характера уравнений и применения факторизации стабилизирующего оператора. В итоге дискретизация Рунге — Кутты правой части для осесимметричных течений привела к дополнительному ограничению шага по времени, причем другая запись исходных уравнений (которая получается в результате дифференцирования левой части уравнения (1) и переноса в правую часть недифференциальных членов) работает при еще меньших временных шагах. Следует отметить, что по соображениям точности для разрешения быстрых процессов в течениях с разрывами приходится
применять эмпирическое ограничение К < 1.1 — 1.3. Таким образом, в расчетах автоколебательных течений дополнительное ограничение не является существенным. Однако оно может оказаться таковым для течений с погранслоями или другими зонами с квазистационарными процессами и мелкой сеткой.
Как указано выше, используемый метод программно реализован для произвольных криволинейных координатных преобразований х = х(<^,п), У = у(<^,п), имеющих необходимую гладкость и отображающих единичный квадрат на плоскости переменных п в криволинейный четырехугольник на плоскости физических переменных х, у. В рамках такого подхода получить хорошие сетки для сложных физических областей бывает затруднительно. Поэтому создана также версия основной программы для случая, когда функции х = х(£, п), У = у(£, п) отображают квадрат с вырезами (£ < <^0,П < По}, (£ > £ъП < П1} на криволинейный четырехугольник с двумя криволинейными вырезами. В описываемых ниже задачах полагалось <^0 = 0, п0 = 0.
2. Обтекание цилиндра с передним выступом
Согласно экспериментальным данным [1], обтекание цилиндра с передним выступом будет нестационарным, если соотношение длины выступа и радиуса цилиндра примерно равно 1. Здесь рассматриваются автоколебательные режимы для случая цилиндрического выступа. Число Маха невозмущенного потока равно 3.7. Используются сетки 91 х 88 и 128 х 124, образуемые расстановкой сеточных узлов на прямолинейных отрезках, соединяющих точки на внешней границе и на теле. Фрагмент расчетной сетки схематически изображен на рис. 1. Постановка граничных условий традиционна: на наветренной части границы параметры потока невозмущены (т. е. головная ударная волна расположена внутри счетной области и не выделяется), на оси симметрии (у = 0) используется соотношение V = 0, на теле — условия непротекания, на подветренных участках свободной границы — "мягкие" условия.
На рис. 2 приводится динамика давления в точке А, находящейся на стыке двух цилиндров (см. рис. 1). Геометрические размеры — Я = 3.2, г = 1, Ь = 3.2 — соответствуют автоколебательному режиму. Расчетные данные приведены для сетки 128 х 124. Можно констатировать наличие почти периодической динамики с периодом 4.03 (на рис. 2, а пред-
Рис. 1. Обтекание цилиндра с выступом. Расчетная область.
ставлены 8 циклов). В расчете на сетке 91 х 88 получен период 3.69. На рис. 2, б приведена часть этих данных в более крупном масштабе. Экстремальные точки отмечены на этом рисунке цифрами 1-6. Структура течения для этих шести моментов времени изображена на рис. 3, 4 в виде распределений линий тока.
По этим рисункам можно сделать вывод, что в потоке происходит периодическое изменение структуры циркуляционных зон. Отметим, в [2], где эта задача решается методом Годунова на расчетной сетке 60 х 70, получено, что в ходе цикла в определенный момент времени отрывная зона разделяется на две, как и в наших расчетах, хотя у нас передняя циркуляционная зона значительно больше. В [2] определен характер процесса как расходов б
Рис. 2. Динамика давления в точке в разных масштабах времени.
Рис. 3. Линии тока в моменты времени 1-3.
Рис. 4. Линии тока в моменты времени 4-6.
ный, что подразумевает периодическое накопление массы газа в циркуляционной зоне с дальнейшим сносом этой массы как единого целого. Наши данные не подтверждают такую картину. Следует также отметить, что в соответствии с другим характером течения в [2] получен и другой период — около 12 (после пересчета к нашей системе обезразмеривания), что примерно троекратно превышает указанное выше значение.
Рассчитано также обтекание конфигурации из двух цилиндров с Ь = 2, остальные размеры и параметры потока не изменились, используется сетка из 91 х 88 узлов. Как и в описанном выше варианте, в этом случае получено нестационарное, близкое к периодическому течение с периодом 3.73 и небольшой амплитудой колебаний давления в точке А, примерно равной 0.2, т. е. почти в 50 раз меньше амплитуды колебаний в случае основной конфигурации. Структура линий тока в этом случае по топологии аналогична изображенной на рис. 3, в и здесь не приводится. Произведен расчет также с Ь = 1.6, при этом получено стационарное течение, снова в целом аналогичное изображенному на рис. 3, в.
3. Натекание струи на плоскую преграду
Из экспериментальных (например, [5]) и расчетных [6] исследований известно, что при натекании струи на плоскую (в частности, безграничную) преграду существуют диапазоны геометрических и струйных параметров, в которых течение носит нестационарный характер. Некоторые из таких автоколебательных осесимметричных течений рассмотрены здесь. Схема течения и расчетная область представлены на рис. 5. Используется сетка, содержащая 181 х 176 узлов.
При постановке граничных условий предполагается, что температура торможения струи равна температуре среды (воздуха с показателем адиабаты 1.4), в которую истекает струя. В итоге на срезе сопла (т. е. на участке от оси симметрии у = 0 до точки с) параметры потока находятся из решения одномерной задачи о течении от точечного источника с заданной температурой торможения и заданными давлением и числом Маха в точке с. Решение этой задачи находится итерационным решением системы алгебраических уравнений в каждой точке данного участка границы. На оси симметрии у = 0 используется соотношение V = 0, на теле — условия непротекания, на свободных границах "мягкие" условия, причем на части верхней границы, прилегающей к стенке, они дополняются со-
Рис. 5. Натекание струи на преграду. Расчет- Рис. 6. Натекание струи на преграду (Ь = 8гс). ная область. Изолинии давления.
отношением Р = 1, в противном случае получалось численное решение с нефизическим падением давления вне струи.
На рис. 6 приведены результаты моделирования взаимодействия струи, истекающей из сопла, с цилиндрическим препятствием. Сопло с углом полураствора а =10° отстоит от препятствия радиуса К = 6.5гс на расстояние к = 8гс. Струя имеет число Маха на срезе сопла Мс = 2, нерасчетность п = 3, показатель адиабаты к = 1.4. Результат моделирования представлен изолиниями давления в развитом течении. Давление в истекающей струе падает и становится значительно ниже давления в окружающей среде, далее, за прямым скачком (отмечен цифрой 1) оно возрастает до значений, которые примерно двукратно превышают давление на срезе сопла. По ходу течения имеется еще один прямой скачок (в некоторые моменты времени он может отсутствовать), давление за которым близко к давлению за скачком, отмеченным цифрой 2.
Проведен также расчет варианта, который отличается от вышеописанного лишь другим расстоянием к = 7.3гс между препятствием и выходным сечением сопла. На рис. 7 показано развитие струи при мгновенном открытии заслонки в выходном сечении сопла. Приводится распределение изолиний давления в моменты времени £ =2.5 и £ = 5. На рис. 7, а за расширяющейся струей имеется область сжатого газа, отделенная с обеих сторон от областей низкого давления скачками уплотнения. На рис. 7, б течение соответствует последовательному сжатию газа в трех скачках уплотнения. Структура давления в развитом течении для этого варианта в целом аналогична изображенной на рис. 6 и здесь не приводится.
Наконец, рассмотрено взаимодействие осесимметричной струи с бесконечным препятствием. Угол полураствора сопла а = 4°, число Маха на срезе сопла Мс = 2.098. Нерасчетность струи равна п = 4.785, показатель адиабаты у = 1.4. Бесконечная стенка находится на расстоянии к = 6.95гс от выходного сечения сопла ортогонально оси струи. На рис. 8 приведена динамика давления в точке торможения. Структура изолиний давления для этого варианта в целом аналогична структуре, изображенной на рис. 6, и не приводится. Однако динамика давления в точке торможения носит более упорядоченный, периодиче-
Рис. 7. Натекание струи на преграду (X = 7.3гс). Изолинии давления в моменты * = 2.5 (а) и * = 5 (б).
Рис. 8. Динамика давления в точке торможения при натекании струи на бесконечную преграду.
ский характер. Можно оценить количество колебаний на рис. 8 примерно как 7, т. е. длина периода т равна 2.29. Если принять в качестве параметров обезразмеривания величины рж = 98 066кг/м/с2, рж = 1.29кг/м3, соответствующие нормальным условиям среды, и радиус сопла тс = 0.015м [6], то частота v, вычисленная по формуле v = (рж/рж)1/2/Гс/Т, равна 8.041 кГц. В [6] приводятся основные тоны колебаний, полученные в эксперименте (8.301, 16.602, 23.715) и в результате численного моделирования (8.921, 17.843, 26.240), т.е. рассчитанная выше частота соответствует первой частоте в этих списках.
Коснемся вопроса об интенсивности пульсаций давления в точке торможения на препятствии. Авторы [6] полагают, что они получили завышенные значения интенсивности вследствие недостаточной диссипативности метода (метода Годунова с линейными распределениями параметров и использованием ограничителей при вычислении коэффициентов этих распределений). Здесь применяется другой метод, однако интенсивность пульсаций для бесконечной преграды такая же, а для преград конечных размеров выше. Поэтому, если завышение интенсивности пульсаций имеет место, то механизм этого завышения, по нашему мнению, может обусловливаться не методом, а осесимметрической моделью течения. Пульсационные волны сжатия при перемещении к оси струи усиливаются ввиду сингулярного характера задачи. В реальном же течении энергия струи распределяется на всевозможные пульсации, не только осесимметрические, и эффект усиления возле оси отсутствует для волн сжатия, отличных от осесимметрических.
4. Взаимодействие солнечного ветра с однородным внешним потоком
Как известно, многие проблемы астрофизики могут исследоваться на основе модели идеального газа. Следует отметить, что при изучении астрофизических задач весьма важно определить характер течения — стационарное оно или нет. Многие численные методы могут не давать достоверной картины: либо подавляют слабую нестационарность вследствие существенной схемной вязкости, либо, наоборот, ее порождают, если генерация энергии в схеме по разным причинам, например вследствие нелинейных эффектов, превосходит ее диссипацию. Приведенный здесь метод имеет высокий порядок, т. е. малую схемную
вязкость, хорошие свойства сходимости по времени, поэтому он может оказаться эффективным инструментом исследования задач астрофизики.
В данном разделе с целью изучения применимости построенной схемы четвертого порядка для исследования астрофизических задач мы рассмотрим достаточно хорошо исследованное взаимодействие солнечного ветра, т. е. сферически симметричного гиперзвукового течения от точечного источника, с межзвездным сверхзвуковым потоком, который будем считать однородным и постоянным.
Сразу отметим, что для исследования астрофизических проблем пришлось использовать новые параметры, управляющие шаблоном схемы. Оказалось, что при больших числах Маха использование формул (7) не позволяет добиться хорошей сходимости. Сделана попытка применения в качестве новых индикаторов "гладкости" разностей давления второго порядка (точнее, их отношения к самому давлению), управляющих адаптивной искусственной вязкостью [11]. Однако в нашем случае оказалось необходимым добавить аналогичный член с плотностью для включения диффузии на контактных разрывах:
d = (1 - )/Pik + ^(6-6+pik)/pik.
С помощью этих параметров вычисляются индикаторы
Yik = min((1 - x)/(1 - dl/d0), 1), x = min(| d | /d0,1), (8)
где d0, dl имеют смысл значений параметра d, при которых индикатор равен нулю и единице. В тестовых расчетах выбрано dl = 0.06, d0 = 0.40. Весовой параметр ф во многих задачах может быть положен равным нулю. Однако в данном случае наличие в потоке контактного разрыва с большим скачком плотности заставляет использовать ненулевое значение ф, так как оно обусловливает включение искусственной диффузии первого порядка на разрыве, которое "размазывает" и стабилизирует его.
Для иллюстрации сходимости по времени схемы с новыми индикаторами гладкости численно исследуется задача об отражении скачка от плоской поверхности для чисел Маха M = 5,10,15, 25. Расчетная область [0 : 2.2] х [0:1] содержит 111 х 45 точек. Скачок падает на плоскую поверхность (описываемую уравнением y = 0) в точке (x = 1, y = 0) под углом в = 35°. Параметры потока перед скачком (справа от линии y = (x - 1) tg(0)) вычисляются по формулам
P =1, р =1, u = -M(k)1/2/sin(e), v = 0. Параметры потока за скачком вычисляются по формулам
P = (2kM2 + 1 - к)/(1 + к), р =((к +1)P + к - 1))/(1 + к +(к - 1)P),
u = M(k)1/2(1 - 1/p)sin(e) - M(k)1/2/sin(e), v = M(k)1/2(1 - 1/p)cos(e).
Падающий скачок, определенный этими формулами, инициирует появление отраженного скачка, постепенно занимающего стационарное положение.
В таблице приводятся данные о числе итераций, необходимых для достижения машинного нуля невязкой R = logl0[maxik | Pik1/рП - 1 I /т\ при различных числах Маха и Куранта. Расчеты проведены при ф = 0.75. Если положить ф = 0, то необходимое число итераций практически не меняется при K = 1 и K =2 и уменьшается на 10-25 % при K =4. Если используются формулы (7) для индикаторов "гладкости", то при M = 10,15, 25 сходимость отсутствует. Вариант M = 25, K =4 не был рассчитан вследствие появления отрицательных значений давления или плотности.
Число итераций, необходимое для установления решения
Число Число Маха
Куранта 5 10 15 25
1 810 730 700 700
2 450 410 430 410
4 1120 1070 1110 -
Взаимодействие солнечного ветра и внешнего потока межзвездной среды исследуется на основе осесимметричной одноатомной модели (к=5/3) [12]. Расчетная область и типичное распределение изолиний давления в потоке приведены на рис. 9. В качестве единицы длины выбран радиус орбиты Земли. Внутренней границей является полуокружность радиуса 15, т.е. она расположена в 15 раз дальше Земли от Солнца. Поскольку параметры солнечного ветра задаются на расстоянии орбиты Земли, для их определения на внутренней границе производится решение алгебраических уравнений одномерной газовой динамики в предположении стационарности и сферической симметричности течения в окрестности Солнца. На других границах используются условия: экстраполяционные соотношения и равенство V = 0 на оси симметрии у = 0, экстраполяционные соотношения справа на подветренной части внешней границы, на наветренной ее части параметры потока фиксированы и равны параметрам межзвездной среды.
На рис. 9 приведены изолинии давления в развитом течении при числе Маха внешнего потока = 2, параметрах солнечного ветра Ме = 5, ре = 10, Те = 100 на расстоянии орбиты Земли, на рис. 10, а — изолинии плотности в этом же течении. Используется сетка 141 х 143. Имеются слабая внешняя ударная волна (линии уровня слева у границы области) и внутренняя ударная волна с большим числом Маха и большим относительным изменением давления. Между двумя скачками имеется контактный разрыв. Его можно увидеть на рис. 10, а и б, содержащих распределения изолиний плотности в моменты времени, интервал между которыми составляет ДЬ = 67.5. Конфигурация внутренней ударной волны на этих рисунках различается в части, максимально удаленной от оси симметрии. Небольшие изменения формы этой части ударной волны продолжаются в течение длительного времени, и стационарной картины течения в расчетах не получается. Слабая нестационарность течения и смыкание контактного разрыва и ударной волны вблизи оси симметрии отличают описываемые результаты от приведенных в [12].
Рис. 9. Взаимодействие солнечного ветра и внешнего потока межзвездной среды. Изолинии давления.
Рис. 10. Взаимодействие солнечного ветра и внешнего потока межзвездной среды. Изолинии плотности.
Эти данные получены при ф = 0.75. Если в формулах (8) использовать значение ф < 0.125, то вначале контактный разрыв становится тоньше, но затем ввиду его неустойчивости на нем появляются волны возмущений (например, на рис. 10, а будут содержаться около двух волн возмущений), в результате чего картина становится существенно нестационарной. При промежуточных значениях параметра £ результаты зависят как от этого параметра, так и от характера начальных данных.
Таким образом, даже базовый этап в решении многих астрофизических проблем — определение структуры взаимодействия солнечного ветра и внешнего потока межзвезной среды — является сам по себе достаточно сложной проблемой вычислительной аэродинамики и требует применения высокоточных и надежных численных методов.
Итак, численные результаты, полученные при использовании неявной матричной схемы четвертого порядка [9], демонстрируют как приемлемое разрешение структуры сложных течений в их динамике, так и отделение стационарных и нестационарных режимов. Вместе с тем в представленных результатах есть некоторые отличия от опубликованных в литературе по данной тематике. Представляется, что в исследуемой сложной области аэродинамики не закончен процесс уточнения структуры основных известных нестационарных течений. Этот процесс идет за счет как повышения точности численных методов и увеличения числа узлов сетки, что особенно важно в расчетах течений с отрывами, так и совершенствования постановок решаемых задач.
Список литературы
[1] Базыма Л.А. Аэродинамика тел вращения с передней срывной зоной. Харьков, 1988. Деп. в ВИНИТИ 3.10.88, № 7262-В88.
[2] Базыма Л.А. О влиянии вдува на стабилизацию течения при обтекании осесимметричного тела с полостью, обращенной навстречу сверхзвуковому набегающему потоку // ПМТФ. 1998. Т. 39, № 4. С. 84-90.
[3] Антонов А.Н., Елизарова Т.Г., Павлов А.Н., Четверушкин А.Н. Математическое моделирование колебательных режимов при обтекании тела с иглой // Мат. моделирование. 1989. Т. 1, № 1. С. 13-23.
[4] Shang J. S., Hankey W. L., Svith R. E. Flow-oscillations of spike-tipped bodies / AIAA. 1980. Paper N 62.
[5] Гапоненко Ю.А., Адрианов А.Л., Безруков А.А., Фаворский В.С. Расчетно-экспериментальное исследование взаимодействия сверхзвуковой струи газа с различными
преградами // Мат. моделирование в механике: Сб. науч. тр. / ВЦ СО РАН. Красноярск, 1997. С. 108-124. Деп. в ВИНИТИ 7.11.97, № 3357-В97.
[6] Адрианов А.Л., Безруков А.А., Гапоненко Ю.А. Численное исследование взаимодействия сверхзвуковой струи газа с плоской преградой // ПМТФ. 2000. Т. 41, № 4. С. 106-111.
[7] Глазнев В.Н., Коробейников Ю.Г. Эффект Гартмана. Область существования и частоты колебаний // ПМТФ. 2001. Т. 42, № 4. C. 62-67.
[8] Пинчуков В.И. О неявных абсолютно устойчивых схемах Рунге — Кутты четвертого порядка // Вычисл. технологии. 2002. Т. 7, № 1. C. 96-105.
[9] Пинчуков В.И. Эффективность неявных схем Рунге — Кутты четвертого порядка в задачах газовой динамики // Вычисл. технологии. 2003. Т. 8, № 6. С. 68-77.
[10] Пинчуков В.И. Трех- и четырехшаговые неявные абсолютно устойчивые схемы Рунге — Кутты четвертого порядка // Журн. вычисл. математики и мат. физики. (В печати)
[11] Jameson A., Schmidt W., Turcel E. Numerical solution of the Euler equations by finite-volume method using Runge — Kutta time stepping schemes // AIAA Paper. 1981. P. 81-1259.
[12] ПогорЕлов Н.В. Исследование высокоскоростных газодинамических и МГД течений: Дис. на соиск. уч. ст. докт. физ.-мат. наук. 2001. 300 с.
Поступила в редакцию 26 августа 2004 г., в переработанном виде —12 ноября 2004 г.
