Anderson-et-al-1 (1185923), страница 58
Текст из файла (страница 58)
Записанные в полярной системе координат уравнения, включая зависящие от времени члены, имеют вид — + — + — +Н=О, дЕ дГ дп дт дз дР (6.34) 11 Д. Азаергон и ар тон ! и последняя точка сетки видится из вершины клина под углом 29.52'. Это должно гарантировать отсутствие интерференции ударной волны при граничных условиях, задаваемых на т1 ., Если используется схема Мак-Кормака, то шаг предиктор можно делать непосредственно от стенки, так как используются разности вперед. Но шаг корректор должен быть модифицирован. Один из способов удовлетворения условию скольжения вдоль поверхности клина состоит в использовании на шаге корректор также передних разностей и последующей записи значения о на стенке с граничным условием о =О.
Несмотря на то что использование разностей- вперед и в предикторе, и в корректоре приводит обычно к неустойчивостям схемы, граничное условие на стенке меняет ситуацию так, что устойчивость решения обеспечивается. Типичные распределения давления при обтекании клина, полученные методом сквозного счета, изображены на рис. 6.9. Эти результаты отлично согласуются с аналитическим решением при т = 1.О, обнаруживая резко выраженную ударную волну, несколько размазанную осцилляциями.
Однако результаты тех же расчетов при числе Куранта 0.7 выявляют дисперсионное поведение, присущее схемам второго порядка, которые обсуждались в гл. 4. 322 Гл. 6. Численные методы решения уравнений течения при этом начало системы координат расположено в вершине клина, а векторы записаны в полярных координатах. Если априори предположить, что течение коническое, то можно получить решение в плоскости Я = сонэ!, исключая из рассмотрения про- 28 И 20 й о и 16 Д 12 0 0.06 ОЛО О.!1 0.12 0.13 0.14 р Рис. 6.9. Распределение давления при обтекании клина, полученное методом сквозного счета; М = 2.0, 6 = 7,6; р — безразмерное давление; — аналитическое решение.
изводные в радиальном направлении. В этом случае приходится решать систему ,7! + 68 + Н = О. дн др (6.35) Система является гиперболической относительно времени, и ее можно решать, пока решение не установится по времени, что даст решение задачи стационарного обтекания клина. В некоторых отношениях зависящую от времени систему уравнений . решать легче, например ее процедура расчета много проще.
Как и в примере 6.2, уравнения движения обычно преобра- зуются в вычислительной плоскости. Одно из наиболее часто 3 6.3. Методы сквозного счета 323 применяемых преобразований предложилн Вивьян !У!ч!апд, 1974! и Винокур ['1г'!покпг, 1974]. Это преобразование (см. гл. 5) гарантирует, что система уравнений в строго дивергентной форме после замены независимых переменных может быть записана в прежнем виде. К его недостаткам можно отнести тот факт, что в преобразованных уравнениях в форме Вивьяна в знаменателе консервативных членов всегда возникает якобнан преобразования. Чтобы избежать ошибок при геометрических преобразованиях, необходимо принимать специальные меры при вычислении метрических коэффициентов. Это подробно будет обсуждаться в гл.
10. гие Рнс. 6.10. Преобразованная область. Сложностей, возникающих при использовании простой прямоугольной сетки в примере 6.2, можно было бы избежать, если бы ударная волна рассматривалась как разрыв. На самом деле в большинстве методов сквозного счета ударные волны на границах выделяются, а внутренние ударные волны улавливаются по мере их возникновения. В то время как одни и те жесоображения лежат в основе методики выделения ударных волн как в стационарных задачах обтекания, решаемых маршевым методом, так и в задачах, зависящих от времени, несколько иная схема иногда применяется для расчета давлениявнутриобласти или за ударной волной, когда исходные уравнения записаны , в дивергентной форме.
Рассмотрим систему уравнений в частных производных, записанную в форме (6.26). Пусть мы используем преобразование (х, у) -ь($, т1), $ = х,- П = у/у,(х), (6.36) где у — у,(х) = 0 есть уравнение, задающее поверхность ударной волны. Как показано на рис. 6.10, физическая область пре.
образуется в вычислительную, в которой ударная волна распо. ложена на поверхности т)= 1.0. Дивергентная форма записи 11* 324 Гл. 6. Численные методы решения уравнений течения Тогда уравнение нормали к ней записывается в виде (6.38) П [1+(иу,Их) 11 Нормальная компонента скорости перед ударной волной дается выражением ц „=и, У = ~ — и — +о ). (6.39) 1 / ока [1+ (Ыу,/Ых)а]ЦЯ ч их Откуда для наклона ударной волны получаем ( ца — ця' — = — ц о ° ию ) и„= цт оа (цт цт ) (цт оа ) (6.40) Величина ца „в (6.40) определяется по перепаду давления на ударной волне, задаваемому уравнением (5.209): (6.41) После того как наклон ударной волны рассчитан, в новом положении известны все величины.
Эта процедура повторяется на шагах предиктор и корректор. Затем опять выделяем ударную волну, считая давление за ней . (или какую-нибудь другую ве- уравнений может быть такой, какую предложил Вивьян, или любой другой формой, в которой имеет место сохранение соответствующих потоковых членов. Мы опять пользуемся маршевым методом для отыскания решения во внутренней части ударного слоя. На ударной волне для получения оценки одной из переменных следует применять одностороннее интегрирование.
Будем считать, что в начальный момент времени нам известно все вдоль поверхности задания начальных данных, включая наклон ударной волны. Далее мы последовательно получаем решение во внутренней части ударного слоя, включая точки ударной волны, продвигаясь в направлении маршевой координаты. Кроме того, интегрируется уравнение для наклона ударной волны с[у,/с[х, что дает уточненное положение ударной волны.
Теперь мы вычисляем наклон ударной волны в ее новом положении, а также другие зависимые переменные, кроме давления. Если давление за ударной волной известно, легко найти плотность и обе компоненты скорости из соотношений Гюгонно— Рэнкина. Требуется получить выражение для наклона ударной волны. Запишем уравнение поверхности ударной волны и — у,(х)=0. (6.37) $6.3, Методы сквозного счета (6.42) где !! — вектор консервативных переменных, а Е и Р— вектор- функции от 1!.
Если схема интегрирования задается уравнением (4.58), то значение 1! на следующем временном слое задается выражением и+~ в И ~(д$3)" ! (дС)"+~~ или Последнее дает второй порядок аппроксимации вектора неизвестных !! ~' на следующем времеинбм слое. Эта схема неявная, поскольку производные 1! да и сама величина 11 берутся на следующем временнбм слое, связывая, таким образом, неизвестные в соседних узлах сетки иа следующем времеинбм слое. Для получения линейного уравнения, которое можно разрешить личину) известной. Этот подход предложил Томас и др.
]ТЬотаз е! а1., 1972]. Так как мы рассматриваем методы решения либо зависящих от времени уравнений, либо уравнений для стационарных сверхзвуковых течений невязкого газа, то эти уравнения относятся к типу гиперболических. Обычно гиперболические системы решаются явными методами. Однако величина шага для большинства явных схем ограничена условием Куранта — Фридрихса— Леви.
Для некоторых задач это может привести к неразумно большим временам счета. Разработка и применение полностью неявных схем для решения гиперболических уравнений — сравнительно недавнее явление. Более ранние попытки такого рода были частично явными или итерационными. Недавно были разработаны неитеративные алгоритмы ([1!пг!епш!!т, К11!ееп, 1973; Вг!!еу, Меропа!д, 1973; Веаш, вагш!пд, 1976]).
Преимущество неявных схем заключается в их безусловной устойчивости. Хотя по сравнению с явными методами приходится делать больше вычислений на один шаг по времени, полное время счета может оказаться меньшим. Мы опишем в общих чертах основную схему Бима — Уорминга !Веаш, вагш)пд, 1976] для уравнений, записанных в дивергентной форме, а затем рассмотрим алгоритм расщепления потоков, предложенный Стегером и Уормингом 16!едег, вагш)пд, 1979]. Рассматриваемая система аналогична уравнению (5.192). Для удобства приведем его здесь: д$5 дв дг — + — + — =О, д! дх дк 326 Гл. 6. Численные методы решения уравнений течения относительно 0"+', используется локальное разложение в ряд Тейлора производных Г и б.
Пусть Е"+' — Е" ] [А](0"+' 13") Г"~' = Г" + [В] (13""' — 13") где [А] = дЕ/д13, [В]= дГ/д(3. После подстановки линеаризации (6.44) в (6.43) для определения 13 "+' имеем линейную систему относительно неизвестной 13 "+'. (6.44) 1[3]+ 2 [ л [А]" + л [В]"3~13" '= =1[3]+ 2 [ л [А]" + л [В]") ~ 13" — Ас], л + л 3 '. (6.45) Прямого решения уравнения (6.45) стараются избегать из-за большого количества вычислений в случае многомерных задач.